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Abstract 

In  this  research,  we  consider  the  transient,  or  time-domain,  scattering  problem 
of  an  overfilled  cavity  embedded  in  an  impedance  ground  plane.  This  problem  is  a 
significant  advancement  from  previous  work  where  more  simplified  boundary  condi¬ 
tions  were  used,  which  can  limit  the  number  of  applications.  This  research  supports 
a  wide  range  of  military  applications  such  as  the  study  of  cavity-like  structures  on 
aircraft  and  vehicles.  More  importantly,  this  research  helps  detect  the  biggest  threat 
on  today’s  battlefield:  improvised  explosive  devices. 

An  important  step  in  solving  the  problem  is  introducing  an  artificial  boundary 
condition  on  a  semicircle  enclosing  the  cavity;  this  couples  the  fields  from  the  infinite 
exterior  domain  to  those  fields  inside.  The  problem  is  first  discretized  in  time  using 
the  Newmark  scheme,  and  at  each  time  step,  we  derive  the  variational  formulation 
and  establish  well-posedness  of  the  problem.  This  sets  the  foundation  for  the  finite 
element  method  used  in  the  numerical  analysis.  Using  both  planar  and  overfilled 
cavity  models,  we  provide  numerical  results  through  the  depictions  of  the  electric 
field  and  radar  cross  section  of  the  cavities. 
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ANALYSIS  OF  TRANSIENT  ELECTROMAGNETIC  SCATTERING  FROM  AN 


OVERFILLED  CAVITY 

EMBEDDED  IN  AN  IMPEDANCE  GROUND  PLANE 

I.  Introduction 


1.1  Motivation 

The  study  of  electromagnetic  scattering  of  plane  waves  has  long  been  an  area 
of  great  interest  in  the  scientific  community.  Many  analytical  and  computational 
methods  have  been  explored,  and  the  various  approaches  to  the  study  of  scattering 
problems  are  matched  by  the  wide  variety  of  applications.  These  include  acoustics, 
optics,  and  microwave  technology,  but  it  is  particularly  appropriate  to  the  military 
regarding  radar  applications. 

At  the  essence  of  any  scattering  problem  is  the  concept  of  the  computation  of 
radar  cross  section  (RCS).  RCS  is  defined  as  a  measure  of  power  scattered  in  a  given 
direction  when  a  target  is  illuminated  by  an  incident  wave,  which  helps  reveal  the 
signature  of  a  target  [62],  This  leads  to  efforts  not  only  to  measure  and  predict  RCS, 
but  also  to  enhance  or  minimize  RCS.  More  demands  and  complexities  in  modern 
design  and  integration  of  systems  require  that  a  solid  foundation  exist  for  all  types  of 
problems  involving  RCS  measurements. 

Intense  research  in  this  area  is  motivated  in  part  by  the  fact  that  of  all  the  con¬ 
tributors  to  RCS,  cavities  have  a  significant  impact  on  a  structure’s  overall  signature 
[5] .  Many  scattering  geometries  have  been  studied  in  the  context  of  military  applica¬ 
tions.  A  common  one  is  the  cavity,  and  examples  of  cavity-like  structures  on  aircraft 


1 


and  vehicles  are  engine  inlets,  cracks,  gaps,  and  cavity-backed  antennas  [47].  Planar 
cavities  representing  these  objects  have  been  studied  in  detail,  but  the  overfilled  cav¬ 
ity  is  important  because  it  may  more  accurately  model  surfaces  that  include  defects 
and  perturbations  [111].  For  example,  this  geometry  is  exhibited  in  aperture  and 
conformal  antennas,  which  are  designed  to  operate  in  the  presence  of  a  ground  plane 
[56], 

For  the  military,  though,  the  battlefield  has  evolved  significantly.  For  example, 
in  today’s  conflicts  and  engagements  in  Afghanistan,  precision  is  critical  not  only 
in  identifying  the  enemy  in  challenging  terrain,  but  also  in  expediting  the  decision 
making  process  through  positive  identification  of  threats.  One  of  the  biggest  threats 
today  are  improvised  explosive  devices  (IEDs).  These  devices  are  concealed  or  par¬ 
tially  buried  in  the  ground  and  are  detonated  by  an  electronic  or  pressure- activated 
trigger  (See  Figure  1).  Information  on  IEDs  is  often  gathered  from  drones  and  sen¬ 
sors,  which  have  the  capability  to  detect  ground  that  has  been  dug  up  and  disturbed 
to  plant  IEDs.  In  addition,  jammers  can  be  used  to  interfere  and  stop  electronic 
detonation  of  IEDs  [74],  Consider  the  following  facts: 

The  number  of  IED  attacks  that  killed  or  wounded  coalition  forces  in¬ 
creased  to  60  in  December  (2009)  from  32  in  December  2008.  The  total 
number  of  IEDs,  including  those  that  were  found  before  they  detonated, 
increased  to  8,690  last  year  (in  2009)  from  3,783  in  2008.  [74] 

The  military  has  improved  its  ability  to  detect  and  eradicate  IEDs  in  2010  through 
the  expanding  use  of  foot  patrols,  planes,  drones,  and  balloons.  This  is  evidenced  by 
a  37%  decrease  in  troops  wounded  or  killed  by  IEDs  [106].  However,  there  still 
will  be  a  need  to  locate  these  devices  in  challenging  terrain.  Clearly,  then,  IEDs 
exhibit  the  characteristics  of  an  overfilled  cavity,  and  their  detection  gives  rise  to  a 
scattering  problem.  In  this  context,  it  is  critical  to  the  military  to  understand  the 
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mathematical  foundation  of  this  scattering  problem  that  is  applicable  in  a  wide  range 
of  environments. 


IEDs  are  used  in  a  range  of  scenarios  by  insurgents.  This  graphic 
shows  one  method  employed  against  coalition  forces. 

♦  1:  Hidden  insurgent  with  detonator  watching  the  road 

♦  2:  Coalition  convoy 

♦  3:  IEDs  buried  in  grass  verge  linked  by  ’daisy  chain’  of  wire 

♦  4:  Anti-tank  mine  used  as  IED 


Figure  1.  IED  Depiction  [53] 


Most  of  the  published  research  and  mathematical  development  for  cavity-type 
problems  generally  assumes  the  ground  plane  is  a  perfect  electric  conductor  (PEC). 
However,  for  this  problem  it  makes  sense  that  an  imperfect  conductor  is  more  phys¬ 
ically  realistic  when  the  surface  in  question  is  the  ground,  which  can  be  represented 
with  an  impedance  boundary  condition  (IBC).  The  impedance  at  the  surface  is  a  mea¬ 
surable  quantity  dependent  on  the  material  defined  on  the  surface,  and  one  would 
expect  this  idea  would  give  rise  to  more  complicated  expressions  mathematically  as 
opposed  to  perfectly  conducting  surfaces  [95].  An  IBC  might  be  evident  on  a  painted 
or  coated  surface  on  a  plane  or  vehicle.  It  is  also  used  to  model  a  honeycomb  or 
perforated  material  [16].  Moreover,  it  is  also  used  in  approximating  the  ground  and 
objects  lying  on  the  ground  (See,  for  example,  [100]).  Thus,  an  IBC  is  most  appro- 
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priate  to  model  IED  detection.  It  is  the  goal,  then,  to  provide  a  solid  mathematical 
foundation  for  this  geometry  with  IBCs. 

Furthermore,  it  is  critical  to  note  that  this  study  ties  directly  to  the  Air  Force  Of¬ 
fice  of  Scientific  Research  (AFOSR)  mission  statement  in  the  research  area  of  electro¬ 
magnetics,  which  is  to  “Conduct  research  in  electromagnetics  to  . . .  evaluate  methods 
to  recognize  . . .  and  track  targets  (including  Improvised  Explosive  Devices)  ...”  [2] . 
Also,  it  falls  in  line  with  Air  Force  Research  Laboratory  (AFRL)  Focused  Long  Term 
Challenges  (FLTC),  specifically  FLTC  3,  which  is  “Dominant  Difficult  Surface  Target 
Engagement  and  Defeat,”  which  “is  focused  on  the  ability  to  deliver  selectable  and 
scalable  non-lethal  or  lethal  effects  against  adversaries  and/or  their  support  activities, 
IEDs  ...  in  an  urban  warfare  environment”  [3]. 

1.2  General  Problem  Statement 

Let  12  C  M2  be  the  cross-section  (cavity  interior)  of  a  ^-invariant  cavity  in  the 
infinite  ground  plane.  That  is,  12  has  the  same  definition  for  all  values  of  z.  We 
will  assume  that  the  cavity  fillings,  with  material  of  relative  permittivity  (er  >  1), 
protrude  above  the  ground  plane.  We  denote  S  as  the  cavity  wall  and  T  the  cavity 
aperture  so  that  (912  =  S  U  T.  The  infinite  ground  plane  excluding  the  cavity  opening 
is  denoted  as  rext,  and  the  infinite  homogenous,  isotropic  region  above  the  cavity  as 
U  =  M+\12.  Furthermore,  let  Br  be  a  semicircle  of  radius  R,  centered  at  the  origin  and 
surrounded  by  free  space,  large  enough  to  completely  enclose  the  overfilled  portion  of 
the  cavity.  We  denote  the  region  bounded  by  Br  and  the  cavity  wall  S  as  12#,  so  that 
12 r  consists  of  the  cavity  itself  and  the  homogeneous  part  between  Br  and  T.  Let  Ur 
be  the  homogeneous  region  outside  of  12#;  that  is,  Ur  =  {(r,  6)  :  r  >  R,  0  <  9  <  n}. 
Refer  to  Figure  2  for  the  complete  problem  geometry. 

The  problem  statement  is  given  the  incident  electromagnetic  wave  ( E \  H1)  im- 
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pinging  on  the  overfilled  cavity,  we  wish  to  determine  the  resulting  scattered  fields 
(. ES,HS ).  This  is  a  two-dimensional  problem,  and  due  to  the  uniformity  in  the  pr¬ 
axis,  it  can  be  decomposed  into  two  fundamental  polarizations:  transverse  magnetic 
(TM)  and  transverse  electric  (TE).  Its  solution  then  can  be  expressed  as  a  linear 
combination  of  the  solutions  to  TM  and  TE  problems. 

We  first  consider  the  TM  polarization  problem,  and  the  following  formulation 
is  modeled  after  Van  and  Wood  (See  [105]).  In  this  case,  the  magnetic  field  H  is 
transverse  to  the  2- axis  so  that  E  and  H  are  of  the  form 


E  =  (0,  0,  E~),  H  =  (Hx,  Hy,  0). 


In  this  case,  the  nonzero  component  of  the  total  electric  field  Ez  satisfies  Maxwell’s 
equations  and  the  following  boundary  value  problem  (see  Appendix  A  for  derivation 
of  impedance  boundary  conditions  in  the  time  domain): 


*  .  ^  d2Ez  n 

—  A.E-  +  er  —  0 


(TM)  { 


d  Ez 
dt 


dt2 
ydEz 
/ 1  dn 


Ez  |t=o  — 


dE7 


dt 


t= 0 


in  fill  U  x  (0,  00), 

on  S  U  Text  x  (0,  00), 
Et0  in  fiUW 


where  £r  =  £/£0  is  the  relative  electric  permittivity,  E0  and  Et  0  are  the  given  initial 

conditions  and  rj  =  \J nr/£r  is  the  normalized  intrinsic  impedance  of  the  infinite 

dEz 

ground  plane.  We  define  the  normal  derivative  — —  =  WE7  ■  n,  where  n  is  the 

on 

outward  unit  normal.  We  are  assuming  that  we  have  a  non-dispersive  material  in  the 
dsr 

cavity,  i.e.  —A  =  0,  or  that  the  permittivity  is  not  a  function  of  frequency,  but  could 
ou 

vary  with  respect  to  position. 


We  observe  the  scattered  field  Esz  solves 


6 


0 


in  U  x  (0,  oo), 


AF,,» 

-AE> +  ~a¥ 


d Esz  77  dEsz 
dt  ji  d  n 


,  dE\  77  d  E  . 


+  -- 


dt  [i  dn 


)  on  rext  u  r  x  (0, 00), 


and  also  satisfies  the  radiation  condition  at  infinity 


lim  \fr 


dE'l  1  dEi 
+  -- 


dr  c  dt 


=  0,  t  >  0. 


(1.2.1) 


The  homogeneous  region  U  above  the  protruding  cavity  is  assumed  to  be  air 
and  hence  its  permittivity  is  er  =  1.  In  U ,  the  total  field  can  be  decomposed  as 
E~  =  Elz  +  E'l  where  El  is  the  incident  field,  and  Esz  the  scattered  field. 

1.3  Function  Spaces  and  Requirements 


The  appropriate  and  relevant  theorems  and  propositions  used  and  cited  in  this 
work  will  appear  at  the  appropriate  points  in  the  manuscript.  Nevertheless,  the  goal 
here  is  to  briefly  discuss  the  motivation  for  the  appropriate  functional  spaces  and 
numerical  methods  that  will  be  used.  We  will  use  the  Lebesgue  scalar  product, 


(u,  v) 


uvdT , 


with  ||u||  =  ^  J  \u\~  dT^j  <  00,  where  T  is  a  bounded  subset  of  M2.  We  refer  to 
these  as  square  integrable  functions. 

We  note  that  Hilbert  spaces,  or  complete  inner  product  spaces,  are  most  applicable 
to  real-world  computational  schemes  [107].  We  will  also  need  to  deal  with  “function 
spaces  that  are  larger  than  the  classes  of  continuous  and  continuously  differentiable 
functions”  [10].  The  reason  is  that  when  dealing  with  the  differential  equation  in  a 
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numerical  framework,  finite  element  subspaces  are  piecewise  continuous  at  best.  Thus, 
we  require  a  less  restrictive  way  in  which  to  define  boundedness  through  the  space 
of  square  integrable  functions,  or  L2(f2/?).  As  previously  defined,  L2({lR)  contains 
sets  of  piecewise  continuous  functions  used  in  the  finite  element  method  (FEM).  To 
ensure  boundedness  of  the  differential  operator  over  the  domain  of  square  integrable 
functions,  the  derivatives  must  also  be  square  integrable,  which  leads  to  the  idea  of 
Sobolev  spaces  [107].  Thus,  the  framework  established  here  lends  itself  to  numerical 
methods,  but  more  importantly,  this  function  space  will  be  needed  for  any  weak 
formulation  of  the  boundary  value  problem  [97]. 

We  first  define  the  Sobolev  space  H1^^)  for  the  bounded  domain  QR  with  norm: 


“II  hhsir) 


u||2+||Vu||2) 


1/2 


where  u  G  L2( QR)  [10]. 

This  can  be  expressed  alternatively,  as  it  is  in  [54],  as: 


H\nR)  =  {ue  L2(nR)  \  ||Vu||22  +  \\u\\2L2  <  oo }. 

Therefore,  the  Sobolev  space  H1^^)  is  the  space  of  functions  in  L2(QR)  whose 
derivatives  are  also  in  L2(Qr).  These  derivatives  are  defined  in  a  weak  sense,  and  this 
space  does  include  continuous  piecewise  linear  functions  [35].  Another  perspective 
is  that  in  order  for  the  test  functions  to  be  well-defined,  we  require  that  the  partial 
derivatives  be  defined  globally  in  nature,  and  “tolerant  of  certain  kinds  of  singular¬ 
ities”  [35].  Derivatives  in  the  classical  sense  are  local  in  nature,  so  we  want  a  more 
global  definition  where  we  can  incorporate  singularities. 

We  will  also  need  to  extend  the  idea  of  Sobolev  spaces  to  the  boundary.  This 
will  be  apparent  when  working  with  the  coupling  of  the  problem  on  the  artificial 


boundary,  Br.  We  define  the  trace  of  a  function  as  the  value  of  a  function  on  its 
boundary.  For  our  problem  this  is  a  linear  mapping  from  the  Sobolev  space  on  the 
domain  Qr  to  the  Sobolev  space  on  the  boundary  Br  [54],  Furthermore,  the  trace 
theorem  (see  Theorem  4.2.2)  states  that  the  trace  of  functions  in  H1(Q,r)  as  mapped 
to  H1/2{Br)  is  well-defined  and  bounded  [10].  We  are  particularly  interested  in  the 
spaces  Hl{yLR),  its  trace  H1/2{Br),  and  the  dual  space  of  H112{Br),  H~1^2(Br).  The 
dual  space  is  simply  the  set  of  bounded  linear  functionals  on  the  space  H1^2{ Br ). 

All  of  this  will  help  lay  the  framework  for  numerical  applications.  It  is  impor¬ 
tant  to  note  that  the  methods  we  will  discuss,  the  FEM  and  the  boundary  element 
method  (BEM)  “belong  to  the  most  used  numerical  discretization  methods  for  the 
approximate  solution  of  elliptic  boundary  value  problems”  [97].  We  will  ultimately 
see  that  a  combination  of  these  methods  will  be  required,  and  is  referred  to  in  the 
literature  as  a  hybrid  finite  element  /  boundary  integral  (FE/BI)  method. 

1.4  Time  versus  Frequency  Domain 

Any  scattering  problem  starts  with  Maxwell’s  equations  which  describe  the  prop¬ 
agation  of  electromagnetic  waves,  and  the  resulting  problem  can  be  analyzed  in  either 
the  time  domain  or  frequency  domain.  Quite  often,  scattering  problems  are  posed  in 
the  frequency  domain  due  to  its  convenience  in  a  mathematical  sense,  and  time  har¬ 
monic  behavior  is  assumed.  Field  quantities  that  are  harmonically  oscillating  with  a 
single  frequency  constitute  time-harmonic  behavior  [55].  Frequency  domain  analysis 
is  “ideally  suited  for  scattering  analysis”  from  plane  waves  in  arbitrary  directions  or 
when  the  scattering  source  is  localized  [56].  One  possible  motivation  for  studying  in 
the  frequency  domain  is  that  in  optics,  most  applications  deal  with  narrow  bands  of 
radiation,  making  studying  fixed- frequency  problems  relevant  [9] . 

Nevertheless,  with  increases  in  computing  power  in  the  digital  age,  time  domain 
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methods  have  increased  in  popularity  [88].  It  is  especially  important  to  the  military 
in  many  strategic  areas  involving  short-pulse  communication  and  radar  systems  [31]. 
It  is  also  important  to  understand  and  predict  behavior  in  the  time  domain  because  of 
the  potential  to  simulate  transient  phenomena,  perform  broadband  characterization, 
and  model  nonlinear  devices  [55].  The  time  domain  “is  ideally  suited  for  antenna 
analysis,  where  one  is  often  interested  in  a  solution  over  a  broad  frequency  band  for 
one  or  a  few  excitations”  [56].  It  is  also  better  suited  for  “visual  representations  of 
understanding  held  interactions”  [88].  Therefore,  research  in  the  time  domain  offers 
many  advantages,  yet  will  be  more  mathematically  challenging.  Furthermore,  a  study 
from  a  mathematical  perspective  will  provide  a  new  foundation  in  the  literature,  since 
most  mathematical  approaches  are  from  the  frequency  domain. 
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II.  Related  and  Previous  Work  /  Literature  Review 


2.1  Solving  the  Electromagnetic  Scattering  Problem 

Electromagnetic  scattering  from  cavities  with  varying  geometries  has  been  studied 
extensively  in  computational  electromagnetics  over  the  years.  Anastassiu  did  a  thor¬ 
ough  review  all  methods  used  for  modeling  electromagnetic  scattering  from  cavities, 
with  a  specific  focus  on  engine  inlets  and  open  ducts  from  aircraft.  He  claimed  this 
was  “One  of  the  most  challenging  problems  in  modern  applied  electromagnetics”  [5] . 
This  could  be  due  in  part  to  the  “appearance  of  spurious  modes  caused  by  interior 
resonances”  [105]. 

Anastassiu  divided  the  approaches  into  exact  or  modal  methods  (such  as  Wiener- 
Hopf  and  mode  matching),  high  frequency  and  spectral  methods,  integral  equation 
methods  (boundary  integral  approach,  method  of  moments  (MoM)),  and  differential- 
equation  methods  (FEM,  FE/BI,  finite  difference  time  domain  (FDTD)).  The  major¬ 
ity  of  these  applications,  however,  are  generally  done  on  perfectly  conducting  surfaces, 
which  may  be  sufficient  and  appropriate  for  most  applications  [5].  Nevertheless,  the 
major  numerical  techniques  in  electromagnetics  remain  the  MoM,  FDTD,  and  FEM. 
It  is  worth  noting  that  hybrid  methods,  involving  a  combination  of  integral  equation 
and  differential  equation  methods,  have  gained  popularity,  particulary  if  it  is  possible 
to  use  a  certain  method  where  it  is  most  computationally  efficient  [107].  From  our 
perspective,  the  FEM  and  BEM  will  be  most  important,  as  we  will  seek  a  hybrid 
FE/BI  approach.  It  is  worth  noting  that  the  BEM  was  developed  in  the  1950s,  and 
has  shown  to  be  a  powerful  tool  in  the  study  of  physical  phenomena  in  an  unbounded 
domain,  such  as  scattering  [89]. 

From  a  mathematical  perspective,  differential  equation  methods  and  integral  equa¬ 
tion  methods  interest  us  the  most.  Purely  mathematical  treatments  of  scattering, 
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however,  appear  to  be  more  limited  in  scope  than  the  numerical  and  computational 
electromagnetic  applications.  It  is  clearly  important  to  establish  the  mathematical 
framework  and  well-posedness  (existence  and  uniqueness)  in  a  general  setting  before 
one  can  implement  these  methods  for  specific  settings  and  geometries.  Colton,  Kress, 
and  Cakoni  have  done  extensive  mathematical  work  in  the  area  of  scattering  theory 
over  the  years,  and  have  used  results  to  transition  to  a  growing  area  of  interest,  in¬ 
verse  scattering,  which  is  the  reconstruction  of  the  object  given  a  scattered  field  (see 
[10]  and  [19]).  Monk  provides  a  thorough  treatment  for  the  underlying  mathematics 
of  the  FEM,  and  discusses  inverse  problems  as  well  [75].  It  is  important  to  note 
that  the  study  of  inverse  scattering  problems  in  the  context  of  IEDs  would  also  be 
a  significant  contribution,  as  particular  shapes  could  be  identified  as  IEDs.  Nedelec 
has  also  contributed  significantly  in  the  context  of  integral  equations  and  integral 
representation  of  solutions  [77].  Angell  and  Kirsch  studied  optimization  methods  in 
the  context  of  antennas,  but  still  showed  the  need  of  direct  scattering  solutions  as  a 
foundation,  using  a  problem  with  IBCs  as  an  example  [6].  Chandler- Wilde  has  also 
done  extensive  research  over  the  years,  focusing  on  scattering  with  IBCs  and  rough 
surfaces  (see,  for  example,  [12]). 

However,  all  of  these  authors’  analysis  is  primarily  in  the  frequency  domain.  All 
have  contributed  extensively  in  this  regard,  with  many  common  threads,  such  as  de¬ 
veloping  cases  involving  IBCs  at  the  surface.  All  use  variational  formulation  methods 
to  establish  well-posedness,  which  provides  a  natural  foundation  for  computational 
methods  such  as  FEM.  Yet  the  complications  involving  an  overfilled  cavity  geometry 
are  not  addressed  by  these  authors.  Therefore,  before  we  consider  the  related  work 
in  this  specific  area,  we  first  consider  the  details  of  the  IBC  and  its  development  over 
the  years. 
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2.2  History  of  Impedance  Boundary  Conditions 


The  impedance  concept  was  first  mentioned  in  1938  for  time  harmonic  field  theory 
by  Schelkunoff  [109].  However,  Shchukin  is  first  credited  with  developing  the  IBC 
by  1940.  Leontovich  also  later  derived  surface  impedance  boundary  conditions  in 
1948,  so  they  all  can  be  considered  the  first  major  contributors  to  the  development 
of  the  impedance  boundary  condition  [81].  These  impedance  boundary  conditions 
were  further  developed  by  Senior  in  1960,  and  they  were  almost  exclusively  used  in 
the  frequency  domain  until  the  early  1990s  [66].  It  was  then  that  these  boundary 
conditions  in  the  time  domain  received  more  attention,  particularly  in  the  area  of 
numerical  implementation.  Therefore,  the  IBC  in  the  time  domain  is  still  a  relatively 
recent  and  evolving  area  of  research. 

2.3  Impedance  Boundary  Conditions  in  the  Time  Domain 

Because  the  impedance  boundary  condition  is  inherently  a  frequency-domain  topic, 
its  conversion  and  implementation  in  the  time  domain  is  difficult  and  remains  an  on¬ 
going  area  of  research.  There  are  two  general  approaches  to  this  conversion:  assuming 
no  frequency  dependence  of  the  material  (i.e.  constant  parameter  materials),  or  as¬ 
suming  the  material  is  dispersive  and  frequency-dependent.  Clearly,  using  a  dispersive 
impedance  boundary  condition  enjoys  the  advantage  of  being  applicable  over  a  larger 
frequency  bandwidth,  but  it  is  more  complex  to  implement  [66]. 

Impedance  boundary  conditions  have  been  used  extensively  in  integral  equation 
formulations,  the  FDTD,  and  the  frequency-domain  FEM.  In  Jin  and  Riley’s  book, 
their  formulations  for  incorporating  the  IBCs  are  initially  done  in  the  frequency  do¬ 
main,  and  then  inverse  Fourier  transformed  [56].  In  fact,  most  of  the  computational 
approaches  in  this  area  seek  to  accurately  model  frequency  dependence. 

Lee,  Shin,  and  Kong  modeled  the  impedance  boundary  condition  in  the  time  do- 
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main,  discretized  using  a  centered- difference  method,  and  used  the  resulting  expres¬ 
sions  in  the  FDTD  method.  They  approximate  the  impedance  boundary  condition 
as  a  rational  function  of  frequency,  accounting  for  the  dispersiveness  of  the  condition, 
and  then  transform  into  a  time  domain  equation  [69]. 

Ida  and  Yufrevev  have  done  extensive  research  into  the  area  of  impedance  bound¬ 
ary  conditions  for  transient  scattering  problems.  Most  of  their  research  involves  for¬ 
mulating  integral  equations  for  use  with  methods  with  the  FDTD  or  BEM  methods. 
They  also  investigate  higher  orders  of  approximation  of  impedance  boundary  con¬ 
ditions,  to  include  derivation  of  the  boundary  conditions  in  time  domain.  Most  of 
their  work  revolves  around  implementation  of  the  IBCs  into  the  time  domain  for  some 
numerical  solver  (See,  for  example,  [52], [115], [114], [116], [113]). 

Chen,  Lu,  and  Michielssen  studied  a  coupled  set  of  time  domain  integral  equations 
with  an  impedance  boundary  condition.  According  to  their  work,  the  method  never 
develops  internally  resonant  fields  and  is  validated  by  numerical  examples  [84], 

This  concept  not  only  is  applicable  in  electromagnetics,  but  also  is  a  popular  area 
in  computational  aeroacoustics  (CAA),  in  which  computation  of  sound  propagation 
through  an  engine  inlet  is  important.  Ozyoruk  and  Long  develop  a  time  domain  IBC 
using  a  Z-transform,  and  they  model  it  as  a  rational  function.  They  incorporate  a 
time  discretized  IBC  into  a  Runge-Kutta  scheme  [80].  Bin,  Hussaini,  and  Lee  develop 
a  broadband  time  domain  IBC  [58].  Maloney  and  Smith  ([72])  and  Schutt-Aine  and 
Oh  ([79])  also  researched  implementation  of  the  IBC  into  the  FDTD. 

Therefore,  the  IBC  has  been  studied  quite  extensively  from  a  computational  per¬ 
spective  in  the  time  domain,  but  is  lacking  from  a  purely  mathematical  perspective. 

Our  goal,  then,  is  to  provide  a  foundation  for  this  process,  keeping  in  mind  we  are  in- 

Qu 

terested  more  in  the  effects  of  a  general  Robin  boundary  condition  (i.e.  — — Y  ku  —  f) 

on 

on  the  structure  of  the  problem,  rather  than  the  explicit  conversion  of  the  IBC  into 
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the  time  domain. 


As  a  result,  we  will  start  with  the  idea  that  there  are  two  general  approaches  to 

conversion.  If  we  assume  constant  parameter  materials,  we  end  up  with  a  straight- 

19 

forward  conversion  where  we  can  substitute  —  for  iuj,  which  resembles  a  first-order 

dt 

absorbing  boundary  condition  [55]. 

There  is  also  another  possible  mathematical  justification  of  this  approach.  If  we 
consider  the  expression  for  surface  impedance  as  in  [8]  as: 


V 


where  a  is  the  conductivity  of  the  surface,  we  note  a  dependence  of  r)  on  frequency 
(cj).  However,  if  we  were  dealing  with  a  situation  where  o  was  proportional  to  some 
multiple  of  frequency,  we  could  state  that  the  frequency  dependence  is  suppressed. 
As  a  result,  our  assumption  would  closely  resemble  the  IBC. 

Finally,  we  mention  that  Wang  in  [108]  studied  transient  scattering  with  a  curved 
absorbing  boundary.  The  numerical  examples  depicting  the  scattered  held  in  this 
work  will  be  important  to  us  during  the  validation  phase  because  FEM  will  be  used 
with  the  Newmark  method.  Thus,  we  can  use  this  study  as  a  baseline  to  ensure  our 
results  compare  favorably. 


2.4  Research  on  Overfilled  Cavities  and  Similar  Geometries 

The  geometry  of  a  cavity  closely  resembles  structures  involving  defects,  pertur¬ 
bations,  grooves,  cracks,  or  rough  or  grated  surfaces,  while  an  overfilled  cavity  can 
represent  partially  buried  or  embedded  objects  in  a  plane.  For  example,  Hamid  and 
Hamid  studied  scattering  by  a  partially  buried  dielectric  sphere  in  the  infinite  PEC 
plane.  They  arrive  at  an  analytic  solution  based  on  the  method  of  images,  and  present 
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numerical  results  for  a  variety  of  conditions  [40].  Also,  Wang  and  Li  studied  a  par¬ 
tially  embedded  cylinder  in  a  random  dielectric  rough  surface.  They  used  the  MoM 
technique  and  provided  numerical  results.  Tsaur  and  Chang  investigated  a  dielectric 
cylinder  buried  in  a  shallow  circular  trough  in  the  two-dimensional  ground  plane,  and 
derived  a  series  solution  as  well  as  numerical  results  in  the  frequency  domain,  using 
a  mode-matching  (or  region-matching)  method  for  both  the  TE  and  TM  cases  (see 
[102]  and  [103]). 

The  study  of  impedance  boundary  conditions  with  objects  or  protrusions  from 
the  surface  can  be  traced  back  to  the  late  1960s.  Goshin  considered  the  problem  of 
an  impedance  cylinder  resting  on  an  impedance  plane,  citing  the  effects  of  surface 
waves  and  expresses  an  exact  solution  [37].  Glisson  studied  an  arbitrarily  shaped 
surface  with  impedance  boundary  conditions  [34],  More  recently,  Swearingen  studied 
acoustic  scattering  from  an  impedance  cylinder  placed  normal  to  an  impedance  plane, 
modeling  scattering  from  trees  in  a  forest  [101]. 

Chandler- Wilde  has  done  a  significant  amount  of  mathematical  research  in  the 
area  of  impedance  boundary  value  problems  of  the  Helmholtz  equation  in  the  half¬ 
plane  (see  [13],  [11],  and  [67]),  as  well  as  with  scattering  from  unbounded  rough 
surfaces  [14],  His  work  has  been  restricted  primarily  to  the  frequency  domain.  Lines 
and  Chandler- Wilde,  however,  did  study  an  inverse  scattering  problem  in  the  context 
of  the  time  domain  with  a  rough  scattering  surface  [70]. 

Mathematical  research  involving  overfilled  cavities  in  the  infinite  ground  plane  has 
primarily  been  limited  to  the  frequency  domain  with  both  an  impedance  boundary 
condition  as  well  as  a  PEC  ground  plane.  Duran,  et  ah,  established  existence  and 
uniqueness  of  “outgoing  solutions  for  the  Helmholtz  equation  in  a  locally  perturbed 
half-plane”  with  impedance  boundary  conditions  [29].  They  developed  a  Green’s 
function  solution  for  the  impedance  half-plane.  In  addition,  they  establish  an  expres- 
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sion  for  the  radiation  condition  due  to  surface  waves  from  the  boundary.  Their  study 
of  the  perturbation  of  the  boundary  resembles  the  geometry  of  an  overfilled  cavity, 
as  they  studied  it  in  the  context  of  outdoor  sound  propagation  or  perturbed  water 
waves  in  a  sea  harbor. 

They  also  separate  the  domain  into  two  sub-domains  connected  through  an  artifi¬ 
cial  semicircular  boundary,  coupled  through  a  Dirichlet-to-Neumann  (DtN)  operator 
(See  Figures  3  and  4),  which  maps  the  values  of  the  scattered  field  on  the  artifical 
boundary  to  its  normal  derivative.  They  establish  results  for  both  the  two-dimensional 
and  three-dimensional  cases  in  the  frequency  domain  (see  [29]  and  [30]).  Their  prob¬ 
lem  closely  resembles  our  problem  geometry,  and  we  will  use  that  approach  as  a 
framework,  particularly  for  the  exterior  problem. 


Figure  4.  Duran’s  Interior  Problem  Geometry  [29] 
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Wood  considered  the  scattering  of  a  time-harmonic  plane  wave  by  an  overfilled 
cavity  in  the  two-dimensional  PEC  ground  plane  [111].  She  also  uses  an  artificial 
boundary  condition  on  a  semicircle  enclosing  the  overfilled  portion  of  the  cavity;  this 
boundary  couples  the  exterior  fields  to  those  inside  the  cavity.  She  solves  for  the  fields 
exactly  via  a  separation  of  variables  (SOV)  technique,  and  then  establishes  existence 
and  uniqueness  of  the  interior  problem  through  a  variational  formulation  using  the 
properties  of  the  DtN  operator  on  the  artificial  boundary.  Huang  and  Wood  extend 
this  work  by  implementing  a  finite  element  method  to  numerically  simulate  the  results 
of  this  problem  geometry  [51]. 

Also,  D.  U.  Kui  also  studied  overfilled  cavities  in  his  dissertation  entitled,  “  Nu¬ 
merical  Computation  of  Electromagnetic  Scattering  From  Large  Cavities”  [65].  Kui 
cites  Wood’s  previously  discussed  work  and  methods  in  this  area,  and  studies  the 
case  for  multiple  overfilled  cavities,  where  the  artificial  boundary  enclosing  all  of  the 
overfilled  cavities  takes  the  form  of  a  semiellipse  instead  of  a  semicircle.  Kui  uses 
the  elliptic  coordinate  system  to  prove  uniqueness  and  existence  for  the  TM  and  TE 
cases  using  SOV  and  the  same  techniques  Wood  uses  in  [111]  for  the  semicircular 
artificial  boundary.  Kui  notes  one  advantage  of  this  semielliptical  artificial  boundary 
for  multiple  overfilled  cavities  is  that  it  increases  accuracy  and  efficiency  due  to  the 
smaller  size  of  the  computational  domain  [65]. 

Jin  and  Riley  propose  the  use  of  the  FE/BI  method  to  analyze  an  antenna  struc¬ 
ture  embedded  in  a  ground  plane  that  partly  protrudes  above  the  surface  (refer  to 
Figure  5).  They  introduce  a  “truncation  surface”  and  employ  image  theory  to  for¬ 
mulate  the  exterior  fields  via  integral  equations.  These  boundary  integral  equations 
take  into  account  the  effect  of  the  ground  plane  on  the  radiated  scattered  fields.  This 
also  eliminates  internal  resonance  problems  [56]. 

The  overfilled  cavity  geometry  has  also  been  studied  from  a  mathematical  per- 


18 


Truncation 


Original 

surface  S„  Mo-  £o 


Image 

surface  5im  Mo-  £o 

(b) 


Figure  3.14  Antenna  structure  partially  residing  in  a  cavity  and  protruding  above  a  ground 
plane,  (a)  Original  problem,  (b)  Equivalent  image  problem  for  the  truncation  surface. 


Figure  5.  Jin  and  Riley  Conformal  Antenna  [56] 


spective  in  the  time  domain.  Van  and  Wood  studied  the  overfilled  cavity  in  the  PEC 
ground  plane  in  the  time  domain.  They  first  discretized  the  problem  in  time,  and 
then  used  the  coupling  technique  and  DtN  operator  to  help  establish  existence  and 
uniqueness  of  the  variational  forms.  In  addition,  they  perform  Enite  element  error 
analysis  and  stability  analysis  for  the  time  stepping  scheme  [105].  Huang  and  Wood 
extended  this  work  by  implementing  a  hybrid  finite  element  and  Fourier  transform 
method  to  provide  an  analytical  and  numerical  solution.  Again,  they  solve  for  the 
fields  in  the  exterior  region  analytically,  and  use  the  FEM  to  solve  for  the  interior 
fields.  The  results  are  finally  inverted  back  into  the  time  domain  [49]. 

Huang,  Wood,  and  Havrilla  use  a  hybrid  finite  element-Laplace  transform  method 
to  analyze  the  scattering  problem  of  an  overfilled  cavity  in  the  PEC  ground  plane  in 
the  time  domain.  Their  approach  is  generally  similar  to  Van  and  Wood  in  [105]; 
however,  they  employ  a  domain  decomposition  technique  integrated  into  the  FE/BI 
method.  Here,  an  overlapping  zone  is  created  between  the  interior  and  exterior  do¬ 
mains  where  the  Laplace  transform  is  used  as  an  analytic  link  condition.  They  also 
provide  numerical  validation  as  well  as  an  error  analysis  [50]. 

It  is  also  important  to  note  that  the  previous  results  assuming  a  PEC  ground  plane 
can  also  be  extended  to  a  perfect  magnetic  conductor  (PMC)  ground  plane.  The  PMC 
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case  is  in  essence  the  dual  of  the  PEC  case.  That  is,  for  the  TM  case,  instead  of  a 


homogeneous  Dirichlet  condition  for  the  tangential  component  of  the  electric  held 

dEz 

(E~  =  0),  it  would  be  a  homogeneous  Neumann  condition  l——1  =  0).  Similarly,  for 

on 

the  TE  case,  instead  of  the  Neumann  condition  for  the  tangential  component  of  the 
magnetic  held,  we  would  have  a  Dirichlet  condition  [94],  Even  though  pure  PMC 
materials  are  not  known  to  exist,  these  results  would  be  an  important  application 
to  the  study  of  metamaterials  which  could  closely  model  this  boundary  condition. 


Applications  in  the  study  of  metamaterials  include  smaller  antennas  and  cloaking 
research  [71]. 


20 


III.  Mathematical  Formulation 


3.1  Approach  to  Solution  and  Semidiscrete  Problem 

In  order  to  solve  a  partial  differential  equation  (PDE),  the  primary  analytic  ap¬ 
proaches  we  can  use  are  SOV,  developing  a  Green’s  function,  and  variational  formu¬ 
lation  [39] .  It  has  been  shown  that  the  separation  of  variables  method  is  appropriate 
for  this  problem  geometry  for  a  PEC  or  PMC  surface.  For  an  IBC,  however,  there 
are  more  complications  involved  with  this  method.  Hanson  and  Yakovlev  address  the 
impedance  plane  problem  through  SOV  by  solving  two  one-dimensional  problems, 
and  combining  the  results  to  get  an  expression  for  the  Green’s  function  as  well  as  for 
the  generated  field  [41],  As  mentioned  previously,  Duran  et  al.  and  Chandler- Wilde 
also  obtained  expressions  for  the  Green’s  function  of  an  impedance  half-plane.  In 
addition,  Politis  et  al.  develops  several  representations  of  the  Green’s  function  for  a 
modified  Helmholtz  equation  with  IBCs  on  the  surface  [82],  Therefore,  a  generalized 
Green’s  function  method  will  be  the  most  appropriate  method  to  solve  the  problem. 
That  is,  we  can  derive  a  Green’s  function  for  an  infinite  planar  surface  with  an  IBC, 
and  for  the  exterior  problem,  get  an  integral  expression  for  the  scattered  field  that 
contains  the  Green’s  function.  Finally,  we  will  use  a  variational  formulation  to  solve 
the  interior  problem.  Thus,  having  explored  all  of  these  primary  analytic  methods 
in  order  to  show  existence  and  uniqueness,  we  will  use  both  the  generalized  Green’s 
function  method  as  well  as  variational  formulation  to  obtain  a  solution. 

Furthermore,  these  analytic  methods  will  provide  the  foundation  for  both  the  finite 
and  boundary  element  methods,  implemented  as  a  hybrid  FE/BI  method.  That  is, 
the  variational  formulation  is  a  framework  for  the  FEM,  which  has  “almost  universal 
applicability”  [97].  In  addition,  “the  use  of  BEM  requires  the  explicit  knowledge  of 
a  fundamental  solution”,  or  Green’s  function,  “which  allows  the  transformation  of 
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the  PDE  to  a  boundary  integral  equation  to  be  solved”  [97].  Furthermore,  “BEM 
are  often  used  ...  to  find  solutions  of  boundary  value  problems  in  exterior  unbounded 
domains”  [97].  Clearly,  then,  these  techniques  are  the  most  appropriate  to  solve  our 
problem.  According  to  Costabel  in  [23],  we  will  also  find  that  the  use  of  boundary 
integral  equations  has  two  important  uses  for  us:  “as  a  theoretical  tool  for  proving 
the  existence  of  solutions  and  as  a  practical  tool  for  the  construction  of  solutions.” 
The  integral  representation  formula  is  a  critical  starting  point  for  these  applications. 

To  formulate  and  solve  this  problem,  we  will  first  decompose  the  entire  solution 
domain  to  two  sub-domains  via  an  artificial  semicircle,  Br,  which  entirely  encloses 
the  overfilled  cavity.  These  two  sub-domains  consist  of  the  infinite  upper  half  plane 
over  the  impedance  plane  exterior  to  the  semicircle,  denoted  Ur,  and  the  cavity  plus 
the  interior  region  of  the  semicircle,  denoted  f 1r.  Refer  to  Figure  6  for  a  depiction  of 
the  two  sub-domains. 

Mr 


The  two  regions  are  coupled  over  the  artificial  semicircle  through  the  use  of  this 
DtN  operator,  thus  exploiting  the  electromagnetic  field  continuity  over  the  artificial 
boundary.  The  variational  formulation  used  for  the  interior  problem  will  allow  us  to 
establish  well-posedness.  In  turn,  the  hybrid  approach  will  come  into  play  as  we  use 
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the  solution  obtained  from  the  FEM  on  the  interior  as  the  solution  on  the  artificial 


boundary.  This  will  allow  us  to  get  an  expression  for  the  scattered  held  in  the  exterior 
domain,  as  desired.  This  technique  of  using  a  semicircular  artificial  boundary  with 
a  DtN  mapping  has  been  used  by  Givoli  and  Vigdergauz  in  [33],  Lin  and  Grosh  in 
[112],  and  Zhao  and  Liu  in  [17].  However,  they  assumed  more  simplified  boundary 
conditions  on  the  planar  surface  allowing  a  SOV  approach  involving  explicit  series 
representation. 

It  is  important  to  first  address  the  possibilities  for  incorporating  the  time  depen¬ 
dence  into  the  problem.  We  would  expect  additional  difficulties  due  to  this  added 
dimension,  where  concerns  regarding  stability  of  the  method  used  become  a  factor 
[22],  The  literature  mentions  three  general  approaches  to  transient  problems  when 
working  towards  formulating  boundary  integral  methods:  space-time  integral  equa¬ 
tions,  Laplace-transform  methods,  and  general  time-stepping  methods  [22], 

For  this  problem,  following  Van  and  Wood  in  [105],  we  will  choose  to  use  the 
Newmark  scheme,  an  implicit  time-stepping  method  that  offers  the  advantage  of 
stability.  This  concept  of  applying  the  time  discretization  directly  to  the  PDE  was 
first  suggested  by  Rothe  in  1930  [15].  This  results  in  a  sequence  of  boundary  value 
problems  for  inhomogeneous  elliptic  equations,  where  the  inhomogeneities  are  simply 
the  solutions  at  the  previous  time  steps  [15].  One  advantage  is  that  once  the  procedure 
is  established  for  the  solution  at  one  time  step,  we  can  “march  arbitrarily  far  in  time 
by  simply  repeating  this  procedure”  [22].  It  is  worth  noting  that  the  Newmark  scheme 
is  primarily  used  in  time-domain  studies  of  dynamic  analysis  of  structures,  though  it 
does  have  some  use  in  electromagnetic  scattering  problems. 

Therefore,  we  will  discretize  the  TM  equations  in  time  by  using  the  Newmark 
time-marching  scheme.  When  compared  to  a  forward,  backward,  or  central  difference 
scheme,  some  advantages  of  the  Newmark  method  are  that  it  can  be  made  unconcli- 
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tionally  stable  and  it  has  the  “best  truncation  error”  for  a  choice  of  parameters  [88]. 
This  method  is  a  two-step  finite  difference  method  in  which  there  is  a  prediction  of 
the  answer  followed  by  a  correction  of  the  predicted  value.  At  each  time  step,  we 
construct  a  nonlocal  boundary  condition  on  the  semicircle  Br  to  couple  the  solution 
in  the  infinite  domain  exterior  to  Br  and  the  solution  in  the  bounded  domain  inside 
Br. 

The  Newmark  scheme  is  defined  by  the  following:  Let  N  be  a  positive  integer,  T 
be  the  time  interval,  St  =  T/N  be  the  temporal  step  size,  and  tn+ 1  =  (n  +  1  )St  for 
n  —  0, 1,  2, ...,  N  —  1.  Denote  the  following  as  approximations  at  t  =  tn+ p 


un+l 

u, 

un+l 

du 

~  dt1 

un+l 

d2u 

~  dt 2 

These  approximations  are  related  by 

un+ 1  =  un  +  Stun  +  ^  [2 f3iin+1  +  (1  -  2 f3)un] ,  0  <  n  <  N  -  1,  (3.1.1) 

un+ 1  =  iin  +  St  [yun+1  +  (1  -  7)hn] ,  0  <  n  <  N  -  1,  (3.1.2) 

where  7  and  f3  are  parameters  to  be  determined  to  guarantee  stability  of  the  time¬ 
marching  scheme  [105]. 

We  will  denote  ul  as  the  incident  held  El,  u  the  total  held  EZ1  and  us  the  scattered 
held  Esz.  The  semidiscrete  problem  is  to  hnd  un+1,  n  =  0, 1, ...,  N,  such  that  we  have 
the  following: 
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Prediction 


~n+l  =  un  +  5tiln  +  _  2$)un ,  (3.1.3) 

itn+1  =  iin  +  5t(  1  -  7 )hn,  (3.1.4) 


Solution 


—A  un+1  +  a2er un+1 

=  oi2erun+l 

in  nR, 

yU  +  l 

rj  dun+1 

on  S', 

p  dn 

un+ 1 

—  us,n+ 1  _j_  yi,n+ 1 

on  Br 

Correction 


(3.1.5) 


iln+1  =  a2(un+1 -un+1), 

(3.1.6) 

un+ 1  =  iin+1  +  St^un+1, 

(3.1.7) 

where  a2 


1 

w 


The  IBC  on  Text  and  S  (See  Appendix  A  for  the  derivation  of  these  conditions) 
becomes: 

n  <7?/n+l 

(3.1.8) 


a,ti  =  _ldun+1 


p  dn 

Using  the  correction  factor  described  above  for  un+1  and  un+1 ,  we  can  express  the 
IBC  in  (3.1.8)  for  the  total  held  as: 


5tja2un+1  +  ''  n'\  =  5tja2un+1  -  iin+1. 

p  on 


(3.1.9) 


Therefore,  the  scattered  held  us,n+1  satishes  the  following  exterior  problem: 
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— Ams,ti+1  +  a2us,n+1  =  a2us,n+1  in  Ur , 

<  Us'n+1(R,0 )  =  g(R,e)  on  BR,  (3.1.10) 

77  <97/S,n~^  - 

(5t7Q;2Ms’n+1  H - — -  =  5t'ya2us’n+1  —  td’n+1  +  A'  on  rext, 

/i  on 

where  g  '=  un+1  —  M*’n+1  and 

Jl  =  -<$t7c*V’n+1  -  - +  5t7a2'Ui,n+1  -  ^’n+1, 

/x  on 

and  the  radiation  condition  is  satished: 

(a  s,n+l  i  \ 

— - +  -iis’n+1  =  0.  (3.1.11) 

or  c  I 

We  note  that  in  considering  the  value  of  K  above,  since  the  incident  wave  is 
known,  the  prediction  will  equal  the  actual  value.  Therefore,  K  =  0. 

3.2  Integral  Representation  of  Solution 

In  what  follows,  we  suppress  the  n  +  1  superscript  from  (3.1.10).  We  seek  the 
solution  for  the  nonhomogeneous  modified  Helmholtz  equation: 

—A u(r)  +  a2u(r )  =  f(r)  in  Ur,  (3.2.1) 

where  r  denotes  position,  f(r)  =  a2us,n+1,  subject  to  nonhomogeneous  boundary 
conditions  of  the  form: 


Au{rs )  +  B 


du(rs) 

dn 


h(rs)  on  Text  . 


(3.2.2) 
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Here,  rs  is  on  the  surface,  and  A  and  B  are  constants  defined  as  A  =  Stja2  and 
B  —  — ,  and  h(rs)  =  8t'ya2us,n+1  —  iis,n+1. 

We  require  that  the  associated  Green’s  function  satisfy  (where  r'  denotes  source 
location): 


— AG(r\r')  +  a2G(r\r')  =  S(r  —  r')  in  U,  (3.2.3) 

dG(r\r') 

AG(r\r')  +  B  ’  =  0  on  Text  .  (3.2.4) 

on 

Using  a  generalized  Green’s  function  method,  we  multiply  (3.2.1)  by  G(r\r')  and 
(3.2.3)  by  u{r )  to  get 


— Au(r)G(r\r')  +  a2u{r)G{r\r')  =  f{r)G{r\r')  (3.2.5) 

— AG(r\r')u(r)  +  a2G(r\r')u(r )  =  <5(r  —  r')u(r).  (3.2.6) 


Next  we  subtract  (3.2.5)  from  (3.2.6)  and  integrate  over  the  enclosed  surface  Ur 
(see  Figure  7)  to  get: 


uR 


(  —  AG(r\r')u(r)  +  Au(r)G(r\r'))  ■  ndS  =  //  5(r  —  r')u(r)dS 

'Ur 


'Ur 

=  u(r') 


f{r)G(r\r')dS 

f(r)G(r\r')dS. 

(3.2.7) 


uR 


We  apply  Green’s  second  identity  in  two  dimensions  to  the  left-hand  side  of  (3.2.7) 
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Figure  7.  Exterior  Domain 


to  get: 


/  (G'(r|r/)Vu(r)  —  u(r) VG'(r|r/))  •  ndi  =  u(r')  —  f(r)G(r\r')dS,  or 

Jc  J JuR 

u(r')  —  If  f(r)G(r\r')dS  +  f  (G(r|r/)Vtt(r)  —  tt(r)VG(r|r/))  •  ndi,  (3.2.8) 
J JuR  Jc 


where  C  is  defined  as  the  contour  enclosing  the  surface,  Ur,  in  Figure  7,  and  C  = 
rext  +  Br  +  Toq.  This  is  a  common  development  and  expected  result  in  the  literature 
(See,  for  example,  [8]  or  [55]). 

Thus,  our  solution  depends  on  the  sources  inside  the  enclosed  region  as  well  as  the 
values  of  u  and  Vu  on  the  contour  [27].  The  Green’s  function,  which  will  be  developed 
separately,  depends  on  the  boundary  conditions  and  radiation  condition.  Therefore, 
we  seek  to  simplify  the  expression  on  the  boundary,  attempting  to  incorporate  the 
impedance  boundary  condition. 

B  du(r)  <9(jt(t* |t*/) 

Following  [26],  we  will  add  and  subtract  the  terms  — — - - - -  from  the 


A  dn 


dn 


boundary  integral  expression  in  (3.2.8)  on  rext  to  get: 


'Text 


G(r 


,^d u(r)  ^  ^dG{r\r')  +  B  du(r)  dG(r\r')  B  du(r)  dG(r\r') 


dn 


dn 


A  dn  dn 


A  dn  dn 


B 


f  i 

rext  A 
du(r] 


(AG(r\r>)*P  _  +  Bd “M  dG^ 


dn 


dn 


dn  dn 


dG{r\r') 


'Text 


dn 

1 

A 


dn 


AG{r\r 


di 


fdu(r ) 
dn 


du(r)  dG(r\r')  .  dG(r\r') 
B  A  7 — X  -  Au(r) 


dn  dn 


dn 


B 


du(r)  dG{r\r') 
dn  dn 


di 


r  ext 


1 

A 


AG(r\r')  +  B 


dG(r\r')  \  du{r 


di 


dn 


dn 


Au(r)  +  B- 


du{r)\  dG(r\r') 


dn 


dn 


di. 


Now  we  use  the  values  of  (3.2.2)  and  (3.2.4)  to  get: 


A 


di 


'  J-  ext 


We  recognize  that  at  the  radiation  condition  will  cause  the  corresponding  contour 
integral  to  vanish,  so  (3.2.8)  further  reduces  to: 


u(r')—  II  f(r)G(r\r')dS  +  /  (G(r\r') V«(r)  —  u(r)VG(r\r'))  ■  ndi 

'Ur  JBr 

1  f  *(r)[®^V  (3.2.9) 


A 


Since  r  and  r'  are  arbitrary,  and  because  the  Green’s  function  is  symmetric  (see 
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(3.3.7)),  we  interchange  r  and  r'  to  get: 


u(r  = 


f  f{r')G(r\r')dS'  +  f  {G{r\r,))^y-  -u(r')dG^ jr  ^)d£' 
uR  Jbr  dn'  dn' 


1 

A 


'Text 


dn'  dn' 

nr')[^IPW,  (3.2.10) 


where  the  normal  derivative  is  taken  with  respect  to  the  primed  (or  source)  coordi¬ 
nates. 


3.3  Green’s  Function  Development 


The  integral  representation  of  the  solution  found  previously  depends  on  finding 
the  associated  Green’s  function.  The  geometry  depicted  in  Figure  7  for  the  exterior 
solution  is  an  annular  sector,  which  accounts  for  the  overfilled  portion  of  the  cavity. 
However,  as  in  Duran  et  al.,  we  first  will  develop  the  explicit  Green’s  function  for  an 
impedance  half-plane  and  apply  it  to  the  modified  problem  geometry.  We  will  note 
that  the  Green’s  function  for  an  impedance  plane  has  been  studied  in  great  detail  by 
several  authors  such  as  Politis  in  [82],  Ochmann  and  Brick  in  [78],  and  most  recently 
Hein  in  [42],  In  fact,  Hein’s  dissertation  significantly  expands  the  work  of  Duran  et 
al.  in  [28],  so  this  is  an  ongoing  area  of  research.  Nevertheless,  we  will  still  follow  and 
expand  the  development  of  the  Green’s  function  in  Duran  et  al.  and  apply  it  to  the 
modified  Helmholtz  equation  (See  [29]  and  [28]). 

We  fix  a  source  point  r'  =  (x',  y')  G  7/  =  \  fh  Thus,  the  Green’s  function  must 

solve  the  boundary  value  problem: 


— AG(r\r')  +  a2G{r\r')  =  8(r  —  r')  in  U, 
,dG(r\r' 


AG{r\r')  +  B- 


dy 


=  0 


on  {y  =  0}. 


(3.3.1) 


Following  Duran  in  [28],  since  no  horizontal  variation  exists  in  the  strict  half-plane 
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problem,  we  can  assume  without  loss  of  generality  that  x'  =  0.  In  order  to  reduce  the 
PDE  in  (3.3.1)  to  an  ODE,  we  take  a  Fourier  transform  in  the  horizontal  direction  x , 
and  denote  the  Fourier  variable  £: 


F[G\=G(£,y)  =  /  G(r\r')e~*xdx, 


dG(£,y) 

dx 

&G(Z,y) 
dx 2 


'  —  CO 
PCO 


=  {-it)  /  G(r\r')e-*xdx  =  {-it)F[G\, 


=  {-it){-it)F[G]  =  - ?F[G ], 


F[S{r  -  r')]  =  F[6(x  -  x’)S{y  -  y')}  =  ~7=S{y  -  y), 

V  Ztt 


S2G((,V)  +  d2G((,  y) 


dx2 


dy2 


+  a2G(£,  y)  =  -7=S(y  ~  y')- 

V  Itx 


Thus,  we  obtain  the  equation  in  the  vertical  variable  y: 


^  ~d  ^dy2  —  +  (^2  +  a2)G(£,y)  =  ^=%-2/') for  y  > 0 

AG(t,  y)  +  f/  =  0  for  y  =  0,  or  rext 

dy 


(3.3.2) 


We  will  solve  this  via  scattering  superposition  (See  [41]),  where  we  seek  a  solution 
of  the  form 

G  =  Gp  +  Gh , 


where  Gp  represents  the  principal  Green’s  function  that  satisfies  the  jump  and  conti¬ 
nuity  conditions,  and  Gh  when  added  to  Gp  will  satisfy  the  radiation  and  boundary 
condition  for  our  particular  problem. 
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First,  we  solve  for  Gp,  away  from  y  —  y': 


d  y2 
dy 2 


+  ({2  +  a2)G»«,!,)=0 
-  (£2  +  a2)Gp(£,y)  —  0. 


The  general  solutions  will  be  of  the  form: 


Gp(£,y)  = 


2+«22/  +  (J2e~ y/Z*+c?y  for  y<y> 

+  for  y  >  y>. 


We  require  that  Gp(^,y)  G  L2(— oo,  oo),  further  reducing  to: 


Gp(e,l/)  = 


-+a2y  fGr  y  <  y' 


C4e  £or  y  >  yf 


Enforcing  the  continuity  conditions  at  y  —  y'  yields: 


'  £2+a2y'  _ 


=  C4e~ 


_  /-t  „-2W£2+a2y' 


C 1  =  C4e 


Finally,  applying  the  jump  conditions  : 


dGP((,y)  v~v'+‘ 

S  ~n~  ,  = 

y  y=y'-e 

-C4y/£2  +  a2e-^S2+aHy'+e)  _  Ciy/£2  +  a2ey/e+<*2{y'-e)  = 


-CAe 


2+ory  _ 


+  OL 2 


32 


Substituting  for  C\  yields: 


-C±e~^^2+a2y'  —  (C^e-2^ 't2+a2yl ^eV^2+a2y' 

-2  CAe-^e+°2yl 


1 


Ch  = 


Ci  = 


V7/2  +  a2 
1 

y/?  +  «2 

gV  €2+cx2y' 

2v/e  +  «2’ 
e— \A2-H*V 


a- 


Therefore,  we  can  express  the  principal  Green’s  function  for  the  impedance  half-plane 


as 


f  g-\/?2+«2j 


Gp(e,l/)=  < 


-  2v^T 


-_eVt2+a2y...y  <  y' 


cr 


j\A2+aV 

[2v^  +  ^e 

f  e\/?2+«2(y-y') 


~\/e+oPy  y  >  yt 


=\ 


■y  <  y' 


or 


3-^/?2+a2(^/-^/,) 


—y  >  y 


which  simplifies  to 


l  2\/  £2  +  ct2 

g— \Jl?+°i?\y' -y\ 

Gp^y)  = 


2x7^  + 


cr 


(3.3.3) 


Now,  to  determine  Gh,  we  first  determine  the  homogeneous  solution  of: 


&Gh(£,y) 
dy 2 


-  (52  +  a2)G'‘K,!/)  =  0, 


which  is 

Gh(i e,  y)  =  +  G6e-^+^^ 

and  C$  =  0  due  to  the  radiation  condition.  Therefore,  we  have  our  Green’s  function 
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representation  as: 


G(i,y)=G’’($,y)  +  Gb($,V) 


e-y/¥+^W-y\ 

2  v/^2  +  «2 


+  C&e~^e+a2y. 


Applying  the  IBC  to  G  above  yields: 


AG{£,y)  +  B 


dG{£,y) 


dy 


p-\/t,2+a?\y' -y\ 

A[  - -  +C6e 


=  0  = 
y= o 

-y/ i2+a2y 


+ 


2x/e  +  «; 

(  x  /p  +  a2e-y/e+<x2\y'-y\  -  r- — -  \ 

B I  . _ = - C6^e  +  a2e-^^y) 

V  2y^  +  ^  V  J 


=  0. 


y= o 


At  y  =  0  we  have: 


/e— y/P+a*\v'\  \  /  fp  i  Q;2g-\/?2+«2h,l  _ \ 

+  C6  +ff(V<:  - =0 

V2VF  +  ^  )  v  2v^T^  v  ) 

p-\J  £,2+a2\y'\  ft 2  I  0,2p--\/$2+q2|2/,| 

A — .  +  AC's  +  . _ = - 

"  '  °  2^2  +  a2 


‘2^ 


-BC6Vf+^  =  0 


^  (A  +  y^2  +  «2)  e-v7^^' 

6_  (A-jBv^2T^2)2VFT^2' 


Now  we  are  able  to  write  onr  final  Green’s  function  representation  as  : 


G(Z,y)  =  G^y)  +  Gh(ty) 

=  e-^+^2\y,-y\  _  (A  +  Bs/e T^2)  e-Ve+^v'  r^/^y 
2^2  +  a2  (A  —  By/  f2  +  a2)  2y/£2  + 
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Factoring  out  \  and  multiplying  by  a  scale  of  ^==  yields: 


l  ( e-Vt2+a2\y'-y\  (. a  +  5^2  +  a2)  e-V«2+°v 

^,y  ~  7^ v ~7/WT7r  ~  {A-B7e+7f) 


e-y/e+*v  j  (3.3.4) 


Now  we  take  the  inverse  Fourier  transform  of  the  previous  expression  to  get  the 
desired  spatial  expression: 


G(r\r')  =  — 
An 


1  f°°  e~'\A2+Q,2h,_d  .  , 


i(x'—x 


V  + 

'  (A  +  Ttf)  \  e-Ve^Hy'+y) 


^.Loo  V  (A-BTeTrf)J 


jW-x)Zd£.  (3.3.5) 


Recognizing  that  the  second  term  in  the  spatial  expression  above  may  be  rewritten 
as 


r  /.  2jh/er^_  \ 

47T ./_«,!  +  A-B^eTrf  v/?+o5 


and  using  the  fact  that  (see  [29]) 


i  r°°  P~\/(,2+a2\y'-y\  a 

/  e*(^  ~x)td£  —  -H^iioiR), 


47T 


a- 


where  R  =  \J {pc'  —  x )2  +  {%)'  —  y )2,  as  well  as  the  fact  that 


we  have 


4  r°°  e— \/i2+«2l2/,-2/l  , 


4vr  y_  oo 

Similarly,  from  [28],  we  have, 


cr 


=  — A"0(o:A). 

27T 


l  T00  g-v/£2+^2h/+y)  , 


47T 


cr 


=  Ko(aR*), 

Zn 
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where  R*  =  \J ( x '  —  x )2  +  (y'  +  y)2.  This  results  in  the  Green’s  function  becoming 


G(r\r') 


1 

2tt 


K0(aR) 


—K0(aR* 
2  n 


e~\/£,2+oi2{y'+y ) 


ei{x'-x)id^  (3.3.6) 


It  is  important  to  observe  at  this  point,  that  this  expression  makes  sense  physically, 
in  that  if  the  surface  was  a  PEC,  then  B  —  0,  causing  the  third  term  in  (3.3.6) 
to  vanish.  This  is  consistent  with  an  assumed  Dirichlct  boundary  condition  on  the 
surface,  and  is  confirmed  in  [104],  Furthermore,  we  note  that  the  first  two  terms  in 

(3.3.6)  correspond  to  the  “classical”  wave  behavior,  whereas  the  third  term  in  (3.3.6) 
incorporates  the  surface  wave  behavior  expected  with  an  impedance-type  surface. 

The  only  question  remaining  is  if  it  is  possible  to  reduce  the  integral  expression  in 

(3.3.6) ,  which  represents  the  surface  wave  behavior,  into  a  simpler  expression  using 
residue  analysis.  In  order  to  facilitate  evaluating  this  integral  using  residue  theory, 
the  third  term  in  (3.3.6)  will  be  rewritten  as 


2  B 
47 r 


r  i(x._x)e 

J-oo{A~B^e+^2) 


dt 


If  £  is  treated  as  a  complex  variable,  then  the  integrand  has  branch  points  at  £ 
and  poles  at  £  =  ±J j  —  a2,  and  it  can  be  shown  the  contribution  due 
poles  is  represented  by 


=  ±ia 
to  the 


A  e-^(y'+y)g\x-*'\\J (#)2-a2 


Hein  makes  the  observation  that  this  expression  “describes  the  asymptotic  behav¬ 
ior  of  the  surface  waves,  which  are  linked  to  the  presence  of  the  poles”  [42],  He  adds 
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that  in  order  to  express  the  proper  physical  results  (that  is,  an  outgoing  progressive 
surface  wave,  and  not  a  standing  surface  wave),  the  above  expression  becomes: 


As  for  the  fourth  term  of  the  Green’s  function  regarding  the  contributions  of  the 
branch  cuts,  we  note  that  Politis  et  al.  in  [82],  Duran  et  al.  in  [28],  and  Hein  in  [42] 
all  reconcile  this  term  in  different  ways,  and  it  is  worth  citing  each  result.  For  example, 
Politis  derives  an  explicit  integral  expression  for  the  branch  cuts  and  further  provides 
a  series  expansion  of  this  term  in  polar  coordinates  [82],  He  states  this  expression 
is  “more  convenient  for  numerical  purposes”;  for  our  case  it  would  be  expressed  as 
follows: 


X  r°°  e£.\x-x'\ 

(f)2  -  a2  +  52>< 

(v7?2  -  a2  COS  ( \/C-  -  a2(y'  +  y))  +  4  sin  (\/{2  -  a2(y'  +  i/)))<f£. 

Duran  in  [28],  however,  reconciles  this  remaining  term  in  the  Green’s  function  by 
using  an  inverse  fast  Fourier  transform  (IFFT)  during  the  numerical  implementation 
to  get  the  spatial  result.  Applying  this  to  our  result  yields  the  expression: 

b_  r 

2vr  J_x  ( A  -  +  «2) 

Hein  states  that  this  term  can  be  reduced  to  be  better  suited  for  numerical  inte¬ 
gration  [42],  Applying  his  technique  to  our  case  yields: 

_ e-sD+D  /  Kl(aJ(x' -x)2  +  e)^=  <%■ 

2?r  U  Vl  ^  W(x'-x)2  +  e 
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Therefore,  Hein’s  recent  development  with  this  portion  of  the  Green’s  function  is 
most  suitable  to  apply  for  both  analytical  and  numerical  methods,  so  we  will  consider 
our  final  Green’s  function  expression  in  the  spatial  domain  for  the  problem  to  be  the 
following: 


G(r\r') 


^ K0(aR )  -  ^K0(aR*) 

Z7V  Zn 

r  A  e~~s(yGy) 

B\/(b)2-“2 


a  Aft  \  fy'+y  / _ 

H - e  bG +v)  /  K\{a\J (x1  —  x )2  +  £2) —  d,£. 

2vr  loo  U  Vl  ^  S  W(x>-x)*  +  e  S 


(3.3.7) 


3.4  Steklov-Poincare  Operator  Analysis 


Recall  the  integral  representation  of  the  solution  in  the  exterior  domain  (the  an¬ 
nular  sector  depicted  in  Figure  7)  is 


u{r) 


f(r')G(r\r')dS' 


1 

A 


h(r' 


dG(r\r') 

dy' 


dx’ 


u(r') 


dnr, 


dO'  :r  e  UR. 


(3.4.1) 


At  this  point,  we  wish  to  implement  an  integral  equation  method  along  the  ar¬ 
tificial  boundary,  Br,  to  couple  the  solution  along  the  artificial  boundary.  We  will 
henceforth  name  the  DtN  operator  as  the  Steklov-Poincare  operator.  This  is  because 
we  will  be  using  integral  operator  theory  to  obtain  an  expression  for  this  operator. 
It  is  expressed  in  terms  of  integral  equations  and  not  series  expressions,  and  the 
literature  is  consistent  in  this  regard. 

First,  let  us  define  formally  the  map  Tr  :  H1G{J3r)  — >  77~1//2(Hr),  where  u  G 

H1/’2{Br)  is  a  radiating  solution  to  the  modified  Helmholtz  equation  mapped  to  its 

du 

normal  derivative  —  e  H~1/2(Br).  Throughout  the  literature,  the  idea  of  mapping  a 
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solution  to  its  normal  derivative  has  several  other  names,  such  as  the  capacity  operator 

[77]  or  the  Calderon  operator  [75].  With  a  Robin  boundary  condition  enforced  at  the 

ground  plane  for  this  particular  problem  geometry,  it  does  not  appear  that  an  exact 

series  representation  is  possible  as  it  was  with  the  PEC  or  PMC  case.  Therefore,  using 

the  integral  representation  of  the  solution,  we  can  express  the  mapping  of  the  integral 

representation  of  the  exact  solution  in  the  exterior  domain  to  its  normal  derivative 

as  restricted  to  Br  as  the  Steklov-Poincare  operator  [54], 

Givoli  further  describes  this  type  of  mapping  as  an  exact  nonlocal  absorbing 

boundary  condition.  The  term  “nonlocal”  refers  to  the  fact  that  we  are  using  a 

boundary  integral  operator  over  the  entire  artificial  semicircle  Br  [32],  The  boundary 

conditions  are  nonreflecting  as  well  in  that  we  would  expect  no  spurious  reflections 

[54],  In  fact,  this  type  of  artificial  boundary  condition  was  the  result  of  many  years 

of  trying  to  reduce  or  eliminate  spurious  reflections  [32],  This  technique  of  truncating 

the  domain  and  using  a  boundary  integral  on  the  artificial  boundary  is  referred  to 

as  a  hybrid  FE/BI  method,  and  is  considered  a  “mainstream  method”  for  analysis 

in  the  computational  electromagnetic  community  [107].  In  addition,  Givoli  makes  an 

interesting  observation  worth  mentioning  in  the  context  of  our  research.  He  states 

that  in  the  area  of  acoustics,  “the  DtN  condition  can  be  perceived  as  a  nonlocal 

Qu 

impedance  boundary  condition ,  relating  the  normal  velocity  (  — )  and  the  pressure 

or 

(■ u )  on  the  boundary”  [32],  He  further  states  that  this  mapping  can  be  difficult  to 
express  in  “closed  analytic  form,”  especially  in  cases  of  “special  geometries”  and  the 
“time-dependent  case”  [32], 
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From  (3.4.1),  we  shift  r  onto  the  artificial  boundary  Br,  to  obtain 


Newton  potential 

A 


UR 


f{r')G{r\r')dS'  — — 


dG{r\r') 


ibr 


(  ~  u(r') 


‘  h(r  /  a  / 

rext  &y 

dG{r\r') 


dx' 


d  nr, 


dnr, 


)<W;reBR. 


Single-layer  potential  Double-layer  potential 


(3.4.2) 


At  this  point,  it  is  worth  discussing  how  to  treat  the  inhomogeneities  incorporated 
with  the  Newton  potential  and  the  term  accounting  for  the  effects  of  the  impedance 
plane.  There  are  several  possibilities  on  how  to  treat  these  terms,  which  result  from 
the  time  discretization  of  the  PDE.  In  our  case,  it  is  divided  into  two  portions,  the 
domain  integral  Ur  as  well  as  the  boundary  integral  rext.  While  using  this  technique 
will  require  further  discretization  of  the  domain  and  a  portion  of  the  boundary,  it  is 
only  done  “for  purposes  of  numerical  integration”  [22],  This  computation  will  ulti¬ 
mately  be  considered  as  the  known  forcing  terms  or  right-hand  side  of  the  variational 
formulation.  As  we  will  develop,  the  other  terms  in  (3.4.2)  will  help  constitute  the 

Steklov-Poincare  operator  expression. 

d 

Taking  a  normal  derivative  — —  of  (3.4.2)  along  Br  yields  the  hypersingular  bouncl- 

on 

ary  integral  equation: 


1  du[r) 

2  dnr 


Normal  Derivative  of  Newton  potential 


=  UILf{r'Mrlr')ds'  y  Lh{r,)8ypdx') 

d  dG{r\r' 


r  ( 

d G(r\r')  du(r') 

'Br  V 

dnr  dnrt 

—  u[r 


dnr  dnrt 


+ 

d6' ;  r  e  Br. 


Adjoint  Double-layer  potential  Hypersingular  operator 


(3.4.3) 


We  define  < p{r' ) 


du(r') 

dnr, 


as  in  Hsiao  et  al.  in  [48],  so  all  of  the  boundary  integral 
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operators  for  r  G  BR  are  expressed  as  follows: 


(S<P)(r)  = 
(Du)(r)  = 
(A<p)(r)  = 

( Hu)(r )  = 


G(r\r')(p(r')dO' ,  (Single-layer  potential  operator) 


(3.4.4) 


ibr 


dGfvlr' ) 

u(r') — - - -d6\  (Double-layer  potential  operator)  (3.4.5) 


'Br  '"*'7* 

f  dG(r\r' 


' Br 


dni 


d  n„_ 

(p(r')dd',  (Adjoint  Double-layer  potential  operator) 

(3.4.6) 

(3.4.7) 


Q  dG(r\r') 

u(r')— - - - -dO'.  (Hypersingular  operator) 


<br 


d  nr  dnmt 


The  mapping  properties  are  respectively  defined  as  follows: 


A  :  H~l/2(BR)  Hl'2(BR), 

D  :  HxI2{Br)  -f  Hx'2{Br) , 
A  :  H~1/2(Br)  ->  H~1/2(Br), 
H  :  H^Br)  H-V^Br). 


Because  we  are  dealing  with  an  operator  of  the  form  —A  +  a2  (the  modified 
Helmholtz  operator  in  Equation  (3.2.1)),  Costabel  states  that  this  identical  form  is 
the  “mathematically  simplest  case”  of  a  formally  positive  second  order  elliptic  partial 
differential  operator  [23] .  Also,  if  a  fundamental  solution  exists  for  use  in  the  integral 
representation  expression,  as  we  earlier  demonstrated  with  the  development  of  the 
Green’s  function  for  the  impedance  half-plane,  then  the  mapping  properties  are  well- 
established  for  the  preceding  boundary  integral  operators.  In  particular,  the  operators 
are  bounded  [21]. 
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If  we  further  define  the  Newton  potential  and  the  impedance  plane  term  to  be 


(Nf)(r) 

( Ph)(r ) 


f(r')G(r\r')dS'  ;  r  G  BR, 


Mr 


1 

A 


h(r 


ndG(r\r')  _ 


dy' 


dx  ;  r  G  Br, 


(3.4.8) 

(3.4.9) 


and  note  that  their  mapping  properties  are 


N  :  L2(UR)  — ►  7/1/2(Hr), 
P  :  L2{Ur)  — >•  H1/2(Br), 


then  N  and  P  are  bounded  operators  as  well  (See  Steinbach  in  [97]). 

Therefore,  denoting  /  as  the  identity  operator,  we  can  rewrite  (3.4.2)  in  operator 
notation  as 

\{Iu){r)  =  (, S<p)(r )  -  (Du)(r)  +  (. Nf)(r )  +  (Ph)(r)  ;  r  G  Br  ==► 

OSVOM  =  +  ^)«W  -  [W)M  +  (-P/i) (»•)] ;  r  eBR. 

Assuming  the  Single-layer  potential  operator,  S,  is  linear  and  invertible  (because  it 
is  H~1/2(J3r)-q lliptic  ;  see  Proposition  3.4.1)  yields 

S-\S<e)(r)  =  +  D)u(r)  -  r‘p/)(r)  +  (Ph){r)]  ;r 

p(r)  =  +  D)u(r)  -  S-‘[(iV/)(r)  +  (Pfc)(r)]  ;re8,.  (3.4.10) 

Here,  as  in  [48],  we  define  7r  as  a  Steklov-Poincare  operator  as  follows: 

(7i«)(r)  =  +  D)u(r), 
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where  we  have  Tn  :  Hl/2(BR)  ->  H-V2(BR). 

We  further  define  the  normal  derivatives  of  the  Newton  potential  and  the  impedance 
plane  terms  to  be 

(N’f)(r)  =  ^:(JJ  f(r')G(r\r')dS1']  ;  r  6  Bn,  (3.4.11) 

(™>(-> -|(- X /_  MrO^V)  ;  r  €  Br,  (3.4,2) 

and  observe  that  their  mapping  properties  are 

N'  :  L2(Ur)  — >•  H~1^2(Br), 

P'  :  L2(Ur)  — >•  H~1/2(Br). 

Then  the  hypersingular  boundary  integral  equation  in  (3.4.3)  can  be  cast  as 

^(/V3)(r)  =  (^)(r)  +  (ifu)(r)  +  (N'f)(r)  +  (P'h)(r) 

{\l  ~  AMr)  =  (Hu)(r)  +  (. N'f)(r )  +  {P'h){r) 

{\l  ~  A )<p{r)  -  (I<p)(r)  =  ( Hu)(r )  -  (lip)  +  ( N'f)(r )  +  (P'/r)(r) 

(-^<d)(r)  -  (4p)(r)  =  ( Hu)(r )  -  (lip)  +  ( N'f)(r )  +  ( P'h)(r ) 

¥>(*0  =  +  A)p(r)  +  (Hu)(r)  +  (N’f)(r)  +  (P'h)(r). 
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Substituting  the  expression  for  (p(r)  in  (3.4.10)  into  the  above  equation  yields 


V(r)  =  (\ I+A)  S-1(i/+BMr)-S-1[(Af/)(r)  +  (Pft)(r)]  +  (Hu)(r) 

+  Wf)(r)  +  (P'h)(r)} 
f(r)  =  (\l+A)S~\l-I+D)  +  H  u(r) 

^ - V* - ' 

Symmetric  form  of  Steklov- Poincare  operator 

-  (\l  +  A)S~'[(Nf)(r)  +  (Ph)(r)]  +  [(J V'/)(r)  +  (F'fcKr)]  ; r  e  B„.  (3.4.13) 

The  above  expression  (3.4.13)  gives  an  alternative  symmetric  form  of  the  Steklov- 
Poincare  operator  as  identified  above  [68].  Therefore,  we  have  two  expressions  map¬ 
ping  u(r)  to  <p(r)  for  r  e  Br. 

fir)  =  +  D)u(r)  -  S~ll(Nf)(  r)  +  (Ph)(r)}  ;  r  €  B„,  and 

fir)  =  (1/  +  A)S-\±I  +  D)+H  u(r)  -  (|j  +  j4)5_1[(JV/)(r)  +  (Ph)( r)] 

+  p7)W  +  (n)(r)];re8s. 

from  which  we  can  deduce  that  the  Steklov-Poincare  operator,  Tr,  may  be  expressed 
as 

Tr  =  S-\h  +  D)  —  {l-I  +  A)S-\l-I  +  D)  +  H.  (3.4.14) 

This  is  a  result  that  is  also  established  in  the  literature  (see,  for  example,  Steinbach 
and  Wendland  in  [98]). 

It  is  important  to  note  here  that  we  are  interested  in  yet  another  symmetric 
representation  of  the  Steklov-Poincare  operator,  Tr,  as 

Tr  =  S~\^I  +  D)  =  S~\h  +  D)SS~ 1  =  S-lZS~\  (3.4.15) 
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where  Z  =  (~I+D)S.  According  to  Hsiao  ,  this  alternative  symmetric  representation 
is  important  “when  coupling  boundary  elements  with  a  symmetric  finite  element 
scheme,  one  can  introduce  a  modified  hybrid  discretization  scheme”  [48].  It  also 
avoids  the  use  of  the  hypersingular  operator.  This  certainly  applies  to  our  problem, 
and  we  will  discuss  this  further  when  we  develop  the  finite  element  approximation. 

Before  we  prove  the  main  theorems  in  this  section,  we  present  several  well-known 
propositions  for  boundary  integral  operators.  The  following  propositions  are  estab¬ 
lished  in  several  sources,  but  here  we  cite  Steinbach  and  Wendland  in  [98]  for  their 
development: 

Proposition  3.4.1  The  single-layer  potential  operator,  S,  is  H~l/ '2  (Br)- elliptic,  i.e., 
(■ Sw,w)l2(Br)  >  cf  \\w\\2h-x/2{Br)  ,Vw  G  H~1/2(Br).  (3.4.16) 

In  addition,  according  to  [98],  since  the  single-layer  potential  operator,  S,  is 
bounded,  we  can  use  the  equivalent  norm  of  the  Sobolev  space  H~1^2(Br)  induced 
by  S  : 

IMI s  =  \J(Swiw)L2(t3R)  Vw  G  H~1/2(Br ), 
as  well  as  the  equivalent  norm  of  the  Sobolev  space  H1!2(Br)  induced  by  S'-1  : 

IMIs-i  =  \J (S^v^) L2{Br)  Vu  G  H1/2(Br). 

Proposition  3.4.2  The  hypersingular  operator,  H ,  is  H1 1 2 (Br)- elliptic,  i.e., 

(Hw,w)L2{Br]  >  c«  ||u>||y/1(Bji)  ,Vu>  e  Hl^(BR).  (3.4.17) 

Proposition  3.4.3  The  double-layer  potential  operator  D  can  be  symmetrized  as  fol- 
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lows: 


DS  =  SA  HD  =  AH  SH  =  -I  —  D2  HS  =-I  -  A 2.  (3.4.18) 

4  4 

Proposition  3.4.4  The  inverse  of  the  single-layer  potential  operator  is  bounded,  i.e., 
cl(S~lw,  w)L2(Br)  <  \\w\\2h1/2(Br)  Vw  G  H1/2(BR).  (3.4.19) 

Proposition  3.4.5  The  Steklov-Poincare  operator,  Tr,  is  H1^2{Br)- elliptic,  i.e., 

(Trw,w)L2{Br)  >  cf  \M\2h1/2{Br)  ,Vw  G  H1/2{Br).  (3.4.20) 


Proposition  3.4.6  The  operator  S  is  self-adjoint,  and  there  exists  the  self-adjoint 
square  root  S1^2  satisfying  S  =  5'1/25'1/2  and  =  ||w||s_i. 

The  following  theorem  is  also  proved  in  [98],  and  due  to  its  importance  in  proving 
Theorem  3.4.8,  we  present  it  here  with  an  expansion  of  the  details  and  reasoning: 

Theorem  3.4.7  For  all  w  G  ^^(Br),  the  following  estimate  holds: 


(-J  +  D) 


w 


—  CK  || W |  c-l  i 


s-1 


where  cr  —  -  +  y  -  —  cf1  ■  cf  >  0. 


Proof  : 


(3.4.21) 
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We  first  note  that 


k,+D ) 


w 


s-1 


(S  +  D)w,  (^/  +  D)w)L2{Br)  (Def’n  of  norm) 

=  ((^J  +  +  D)w,w)l2(br)  (Symmetry  in  Eq  ref) 

=  (Trw,  w) l2(br)  ~  ( Hw ,  w)l2{Br).  (3.4.22) 


Now  we  further  establish  the  properties  of  the  inner  products  as  in  [98], 


(Trw,w)l2(br)  —  (S  1^2STrw,S  1^2w)l2(br) 


<  ||s-1/2srRy|MBji) 
=  iistr^iIs-i  IMI  s-i 


||  S~1/2w\ 


L2(Br) 


C2i  +  d) 


w 


s-1 


Fils-1 


(Hw,w)L2(Br)  >c1  \\w\\h1/2(Br) 

>  cfcf(5-1w,w)L2(Bfl) 


=  cici  IMIs- 


Then  we  can  rewrite  (3.4.22)  as 


C,I  +  D) 


w 


s-1 


-  ( Trw,w)l2(Br )  -  (Hw,  w)l2(Br) 


< 


2I  +  D ) 


w 


s-1 


I  H  S  II  II 2 

s_i  —  c1  cx  |Flls-i 


(3.4.23) 
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If  we  let  a  =  (|/  +  D)w |  S-1  >  0  and  b  =  ||w||s_i  >  0,  then  we  have 


1 

2 


a2  <  ab  —  cf  cf b2 
2 


a  a 
V  ~b 


,a\2  a  < 

V  6-  Ciq 


cfcf 

H„S 


& 


JT  S 


6  +  4  -  4  Cl  Cl 

®  _  1\2  <  1 

6  2'  “  4 


Ci  c 


i  'i 


1  IT  a  a  1 

_ r 6  <  _  <  - 

4  1  1  S  f,  S  2 


4  flCl' 


Therefore  we  have 


\I  +  D) 


w\ 


s-1 


\w  C_1 


<  l  +  J\-c?cf  =  cKn 


Before  we  frame  the  variational  problem  in  the  next  chapter,  we  prove  the  following 
theorem,  similar  to  Cakoni  and  Colton’s  Theorem  5.20  in  [10]: 


Theorem  3.4.8  The  Steklov- Poincare  operator,  Tr,  is  a  bounded,  linear  operator 
from,  H1/2{Br)  — >  H~1/2(Br) .  Also,  the  principal  part  of  Tr,  referred  to  as  Tr,p, 
satisfies  the  coercivity  estimate 


(Tr,pU,u)  >  C  11^11^1/2(5^) 

for  some  C  >  0,  such  that  the  difference  Tr  —  Tr,p  is  a  compact  operator  from 
H^iBn)  ^  H-^Br). 

Proof  : 

To  prove  boundedness  of  the  Steklov-Poincare  operator,  Tr,  we  use  the  equivalent 
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norms  of  the  Sobolev  spaces  defined  earlier,  as  well  as  Proposition  3.4.6.,  which  yields 


\\TRw\\s  -  yj  ( STrW,Trw)L2(Br) 

=  \J  (SS-^I  +  D)w,  +  r>MMB„) 

=  \j  +  D)w, 

*#+HL 

=  (-/  +  -D)w  <  c#  ||w||5_i .  (3.4.24) 

2  c — i 


Linearity  follows  directly  from  the  expression  in  (3.4.14). 

At  this  point,  having  shown  boundedness  and  linearity,  we  define  the  principal 
part  of  Tr  to  be  7r,p,  which  corresponds  to  the  Laplacian  portion  of  the  operator. 
Ideally,  for  our  problem,  this  would  be  the  case  where  a2  =  0  in  (3.3.1).  How¬ 
ever,  because  we  have  defined  a2  =  ,  we  note  that  this  does  not  correspond 

(ot)  p 

to  a  physical  situation,  but  is  only  studied  in  order  to  help  frame  the  variational 
formulation  problem.  According  to  Duran  et  al.  in  [29],  the  principal  part  of  this 
operator  can  be  associated  with  the  Green’s  function  for  the  Laplace  equation  with 
a  Neumann  boundary  condition  enforced  at  the  ground  plane  (i.e.,  referring  again  to 
(3.3.1),  a2  =  0  implies  A  =  0  and  that  the  normal  derivative  of  the  Green’s  function 
must  be  equal  to  zero  at  the  ground  plane).  The  Green’s  function,  as  shown  in  Duran 
in  [29],  is 

G(r\r')  =  log-R  +  log-R*. 

Z7T  ZTT 


In  fact,  this  idea  is  supported  by  Khoromskij  and  Wittum  in  [60],  where  they  seek 
a  simpler,  alternate  expression  for  the  Steklov-Poincare  operator  previously  derived 
in  (3.4.14).  They  derive  an  explicit  integral  representation  for  the  approximation  of 
the  Steklov-Poincare  operator  whose  principal  part  is  the  related  Green’s  function. 


49 


Furthermore,  this  explicit  representation  is  in  hypersingular  form.  This  idea  is  also 
explored  in  [57]  and  [59]. 

Duran  goes  on  to  state  the  following: 


That  Green’s  function  has  been  obtained  by  a  principle  of  symmetry  with 
respect  to  the  x\—  axis  (i.e.  x— axis).  Thus,  the  continuity  and  coercivity 
properties  of  its  associated  integral  operator  (i.e.  the  principal  part  of  the 
operator)  on  S+  (i.e.  the  half-circle,  or  onr  semicircle  BR)  are  reduced  to 
the  whole  circle  S.  [29] 

This  makes  sense  in  the  context  of  the  Laplacian  in  that  it  is  radially  symmetric  [73] . 
Thus,  we  can  study  the  properties  of  the  case  where  a2  =  0  for  a  disk,  denoted  £>disk, 
which  is  comparable  to  our  operator  Tr,p  extended  to  the  whole  circle.  We  note  that 
if  u  is  a  radiating  solution  outside  a  disk  £>disk  of  radius  R,  then  we  can  separate 
variables  and  obtain 


OO 

u(r,  6)  =  ^  anr~neind,  r  >  R  and  0  <  9  <  2i r. 

—  OO 


As  in  Cakoni  and  Colton  in  [10],  we  can  state  that  for  data  on  the  disk 


OO 

“lasmsK  =  S^bnein6] 

— OO 


with  coefficients  bn 


anR  n,  the  mapping  of  TR  P 


dii 

:  u  — >  — —  yields 
on 


Tn.pu  —  —  jjbnein0- 


Then  it  follows  as  in  Cakoni  and  Colton  that  for  a  disk  £?disk5 


'd£>DISK 


Tn,Pwwds  >  C  \\w\\2Hl/2(dB 


disk) 
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Then  through  an  extension  argument  we  can  establish 


T, 


'dt3R 


R  pwwds  >  C  ||ic || ■ 


It  can  then  be  shown  from  several  sources  (see,  for  example,  [29]  or  [10])  that  Tr—Tr^p 
is  a  compact  operator  from  H1^2{ Br )  — >•  H^1^2(Br).  These  properties  are  also  more 
explicitly  stated  and  proved  in  Kirsch  in  [61],  as  well  as  in  Hettlich  in  [43].  □ 

As  a  result,  this  theorem  will  help  us  prove  the  well-posedness  in  the  next  chapter. 
Even  though  we  do  not  have  an  explicit  series  representation  for  the  operator,  we  used 
an  important  idea  in  that  the  principal  part  of  the  Steklov-Poincare  operator  satisfies 
a  coercivity  estimate.  This  will  be  discussed  again  during  the  variational  formulation 
chapter. 


3.5  Green’s  Function  Analysis:  Derivative  and  Singularity  Evaluation 


Having  both  the  integral  representation  of  the  solution  in  the  exterior  domain 
in  (3.2.10),  as  well  as  the  Green’s  function  in  (3.3.7),  we  are  ready  to  establish  the 
existence  of  the  normal  derivatives,  which  will  be  required  for  numerical  implementa¬ 
tion.  Singularities  will  also  need  to  be  addressed  during  the  boundary  integral  matrix 
construction. 

The  normal  derivative  of  the  first  two  terms  of  the  Green  function  as  expressed  in 
(3.3.7)  is  given  by: 


c) 

dnr, 


r  1 
.2tt 


K0(aR) 


—K0(aR* 

Ztt 


a 


(x'  —  x,y'  —  y ) 


x)2  +  (y'  -  y) 

°  K1(aR‘)  ,  + 

2vr  a/  Jx'  —  x)2  +  ( y '  +  y)2 


=  ■  V 

•  nrt 


where  rL,/  =  ,  .,,,  which  is  the  outward  unit  normal  vector  to  the  boundary  at  the 
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point  r' .  At  this  point,  we  address  how  to  evaluate  a  singularity  for  the  normal 
derivative  evaluated  at  these  first  two  terms  of  the  Green’s  function.  Specifically,  we 
will  be  interested  in  evaluating  the  normal  derivative  as  R  and  R*  both  approach 
zero.  We  will  derive  the  case  for  R  — >  0  and  R*  can  be  derived  analogously. 

We  first  note  that 


cos  9  = 


(x'  —  x,y'  —  y ) 


Vix1  ~  x)2  +  (y‘  -  y)2 


■  nr> 


where  6  is  the  angle  between  —  n  t  and  r  —  rf,  and  this  is  a  common  vector  calculus 

result  that  can  be  derived  from  the  law  of  cosines.  We  also  use  the  fact  that  for 

1 

small  arguments,  K0(x)  ~  —  loghr),  so  it  follows  that  KRx)  ~  — .  It  is  also  noted  in 

x 

1 

Wiersig  in  [110]  that  cosd  =  -y  |r  —  r'|,  where  y  is  the  curvature  of  the  curve  at  the 
boundary. 

Therefore,  taking  the  limit  as  R  — >  0,  we  have 


lim 

R-^  o 


2n  J  yj  (x'  —  x)2  +  (y  —  y )2 


a 


n™)  V  2ttJ  \aRJ 

=  lim  (-—)(—) 
R-> o  V  27 tJ  \aR/ 

=  lim  — y. 

R^O  47T 


nrt 


2X|r 
\xR 


However,  the  curvature  y  is  finite,  and  furthermore  is  0  for  straight  line  segments, 
which  we  will  be  using  for  the  boundary  element  method,  so  we  can  conclude  that 


lim 

o 


(x'  -  x,  y'  -  y) 

A /<y  -  x)2  +  (y1  -  y)2 


■  hrr  =  0. 


which  we  will  need  to  evaluate  the  expression  at  a  singularity.  As  previously  stated, 
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this  will  also  apply  evaluating  the  case  when  R*  — >  0. 

Also,  the  previous  result,  although  not  explicitly  proven  as  we  have  done,  is  used 
in  such  works  by  Schneider  and  Salon  in  [91]  and  [93]  and  [90],  and  Koopmann  and 
Benner  in  [63].  They  implement  these  results  in  the  boundary  integral  formulation 
phase  in  dealing  with  singularities  of  modified  bessel  functions  of  the  second  kind. 

The  derivatives  of  the  third  term  of  the  Green’s  function  expression  in  (3.3.7)  can 
be  expressed  as  follows: 

d  G3 
dy 

dG3 
dx' 
d  G3 
dx 


The  derivatives  of  the  fourth  term  of  the  Green’s  function  expression  in  (3.3.7) 
can  be  expressed  as  follows  for  y  and  y': 


dGi  _  dG4 
dy  dy' 


B  2tt 


— Kx(a y/ (x'  -  x )2  +  (y'  +  y)2)- 


y'  +  V _ 

\x'  —  x)2  +  (y'  +  y)2 


The  derivatives  of  the  fourth  term  of  the  Green’s  function  expression  in  (3.3.7)  can 
be  expressed  as  follows  for  x  and  x'\ 
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Once  again,  in  dealing  with  the  singularities  for  these  expressions  in  the  boundary 
integral  formulation  phase,  we  can  refer  to  the  works  cited  previously,  or  Ramesh  and 
Lean  in  [86]. 

At  this  point,  we  can  discuss  an  alternate  expression  for  the  fourth  term  of  the 
Green’s  function  for  more  effective  computation.  This  is  done  in  order  to  “isolate  the 
poles”  and  “adequately  treat  the  slow  decrease  at  infinity”  when  y'  +  y  —  0,  which 
would  occur  when  on  the  impedance  plane  (Text)  [42].  For  the  case  when  |^|  >  |a|, 
we  can  apply  Hein’s  technique  to  our  case  and  get: 


where  Ei  denotes  the  exponential  integral  function. 
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It  follows  that  the  derivative  with  respect  to  x'  of  the  preceding  term  is 


and  with  respect  to  x  is 


These  are  obviously  complicated  terms  to  analyze  and  implement,  and  we  leave  the 
computations  with  respect  to  the  y'  and  y  variable  to  the  numerical  phase. 

In  summary,  the  normal  derivatives  of  the  Green’s  function  can  be  evaluated  at 
any  point,  provided  the  gradient  can  be  calculated,  by  observing  that 


dG 

/ dG  dG 

dnr, 

\dxr  dy' 

8Cr 

/ dG  dG 

dnr 

V dx  ’  dy 

•  nr', 

■  hr- 
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Furthermore,  we  mention  that  the  hypersingular  operator  will  require  the  computa¬ 


tion  of 


d  dG(r\r') 


.  This  will  require  the  values  of 


d2G  d2G  d2G 


-,  and 


uxwxx  wx  .  xxxxu  vvxxx  xv^m-xx^  txxo  »wxuou  ux  ry  .  )  ry  O  O  O  O  / 

onr  onr /  oxox  ayox'  oxoy'  oyoy' 

to  be  computed.  These  can  also  be  computed  during  the  computational  phase,  but 


care  is  needed  in  evaluating  the  singularities  and  ensuring  the  proper  expression  is 


obtained. 
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IV.  Variational  Formulation 


4.1  Recasting  the  Boundary  Value  Problem 

Since  the  ultimate  goal  of  the  previous  mathematical  formulation  was  to  lay  the 
foundation  to  establish  existence  and  uniqueness,  in  this  section  we  find  a  weak  so¬ 
lution  of  the  variational  form.  We  first  note  that  on  the  artificial  boundary,  Br ,  the 
normal  derivative  of  the  total  electric  field  satisfies  the  continuity  condition: 


du 

dn 


r=R 


dul 

dn 


r=R 


dus 

dn 


(4-1.1) 


r=R 


Recall  the  expression  derived  previously  in  (3.4.10): 


<e(r)  =  +  D)u(r)  -  S-'[(Nf)(r)  +  (Ph)(r)]-,  r  e  Br, 


where  the  Newton  Potential  and  impedance  plane  terms  are 


(Nf){r)  —  //  f(r')G(r\r')dS'-,  r  e  Br, 

J  JuR 

( Ph)(r)  =  ~f  h(r')9G ^  ^ dx r  e  Br, 


'Text 


dy' 


(4.1.2) 


(4.1.3) 


where  f(r)  =  a2us,n+1  and  h(rs)  =  5t^a2us’n+l  —  us,n+1.  This  allows  us  to  write 


(Vus’n+1)(r)  =  a2  /  /  us,n+1G(r\r')dS' ;  r  e  Br, 

’Ur 


(4.1.4) 


1 


(P«s,n+1)(r)  =  -- 


r  ext  L 


8t^a2us’n+1  -  iis'n+1 


reBn.  (4.1.5) 
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Then  we  can  define 


($Rfis’n+1)(r)  =  ( Nus’n+1)(r )  +  (Pus,n+1)(r), 


and  note  that  since  the  Newton  Potential  and  impedance  plane  terms  are  bounded, 
then  is  bounded.  We  can  also  define 

{^Rus’n+1)(r)  =  S-\^us'n+1)(r), 

and  note  that  since  S'-1  is  bounded  (see  Proposition  3.4.4.),  then  T R  is  bounded. 
Now  (3.4.10)  becomes 

<p{r)  =  =  (Tru){v)  -  (tf Rus,n+1)(r)f,  r  e  BR. 

onr 

Applying  this  to  the  continuity  property  in  (4.1.1)  yields 

+  (Tr(u  -  u1))  -  (yRus’n+1)-,  r  e  Br.  (4.1.6) 

r=R 

This  is  the  transparent  boundary  condition  on  Br,  so  we  may  now  rewrite  (3.1.5)  as 

— Awn+1  +  a2£run+1  =  a2£run+1  in  QR, 

<  5t^a2un+1  +  —  — — -  =  <5t7o;2-un+1  —  iin+1  on  S,  (4.1.7) 

H  on 

fhin+ 1 

__  _  rRun+l  =  —  -  ^Rus'n+1)  -  TniiP)  on  Br. 

Now  we  seek  to  solve  (4.1.7)  through  a  variational  method. 


du  dul 

dn  n  dn 
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4.2  Variational  Formulation 


Instead  of  enforcing  the  boundary  conditions  on  the  test  function  space  V  as  in 
[105],  we  choose  to  define  the  subspace  V  simply  as  H1( £Ir),  which  has  the  i71-norm, 
i.e., 

\\u\\v  =  II^IIh1^)  • 

The  variational  formulation  of  (4.1.7)  will  be  to  find  such  that 

brM(u,  v)  —  F(v)  Vu  G  V.  (4-2.1) 


Note  that  because  we  are  enforcing  Robin  boundary  conditions  on  the  surface  S,  these 
boundary  conditions  are  considered  natural  and  do  not  need  to  be  incorporated  into 
the  solution  space  V.  That  is,  the  boundary  conditions  “appear  naturally”  and  are 
automatically  incorporated  into  the  boundary  integral  expressions  and  are  satisfied 
by  weak  solutions  [46].  We  first  must  construct  the  sesquilinear  form  &tm(«,  v)  as 
well  as  the  conjugate  linear  functional  F(v).  To  do  this,  we  multiply  (4.1.7)  by  a  test 
function  v  G  V  to  obtain  (we  suppress  the  n  +  1  superscript): 


/  A uvdxdy  +  a2  /  eruvdxdy  =  a2  eruvdxdy 

J  Qr  J  Qr 

r]  f  du 


dt'ya 2  /  uvdi  H —  /  7 —vdi  =  dt^a2  /  uvd£  —  /  iivdl 
Js  y  Js  on  Js  Js 

(  Tnuvdi  —  [  ^vdl-  I  yRusm 


^vd£ 


ibr 


dn 


ibr 


>br 


dn 


'  Br 


'Br 


(4.2.2) 

(4.2.3) 
TR{ul)vdl 

(4.2.4) 


Noting  Green’s  identity: 


/  Vit  •  Wvdxdy  =  /  —vdi  —  [ 
'nR  Jbrus  on  JQr 


Auvdxdy, 
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we  can  rewrite  (4.2.2)  as 


Vu  ■  'Vvdxdy 


[  *w~vdl  +  a2  [  eruvdxdy  =  a2  [  eruvdxdy ,  or  (4.2.5) 
IbrUS  On  JnR  JnR 


Vu  •  Vvdxdy  —  (  f  ^-vdi  +  [  ^vdi]  +  a2  [  eruvdxdy 

on  Jq  on  )  ,/On 


^ Or 

=  a2  eruvdxdy. 

J  ftR 


(4.2.6) 


Now  substituting  the  known  values  from  (4.2.3)  and  (4.2.4): 


[  ^ m  =  a 

/s  dn  77 


dt'ya2  /  uvdi  —  /  unfit  —  dt'ya2  /  uvdi 

Js  Js  Js 


,  and 


f  —  vdi  =  /  Tnuvdl  +  f  ——vdi  —  f  4 >  RUsvdi  —  f  Tr{ux)v(M. 

Ibr  on  JBr  JBr  dn  JBr  JBr 


and  letting 


dul 

dn 


-  Tru\ 


r=R 


we  obtain  for  (4.2.2): 


V«  ■  Vvdxdy  — 


in  F 


/  Tfiuvdi  +  /  Jvdi  —  /  'k  Rusvdi 
IBr  Jbr  Jbr 


+  —dt^a"  /  uvdi - /  uvdi  —  —dt'ya2  /  unfit' 

V  Js  V  Js  V  Js 


+  a2  £ruvdxdy  =  a2  £ruvdxdy. 

J  £Ir  J  £Ir 


(4.2.7) 


Observing  that  for  77  ^  0,  gathering  appropriate  terms  yields 


/  Vu  ■  Vvdxdy  —  /  TRUvdi  + —dtja2  I  uvdi  +  az  I  £ruvdxdy 
’nR  Jbr  V  Js  JnR 

=  /  Jvdi  —  /  d' Rusvdi  +  —dt'ya2  /  uvdi  —  —  uvdi  +  a2  £ruvdxdy. 
Jbr  Jbr  V  Js  V  Js  JnR 
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So  now  we  define  the  sesquilinear  form 


i>TM(u,v)  —  /  'Vu-'Vvdxdy—  /  Tpuvdi  + —St^a2  /  uvdi  +  a2  /  eruvdxdy , 

Jbr  V  Js  JnR 

(4.2.8) 

as  well  as  the  bounded  conjugate  linear  functional  term 


d  sr+ _ 2 


F(v)  —  /  Jvdi  —  /  +  —  dtja2  I  uvdi  —  —  I  uvdi  +  ad  /  eruvdxdy , 

Jbr  V  Js  V  Js  JnR 

(4.2.9) 

and  we  follow  [10]  to  rewrite  brM(u,v )  =  &tmi(«4')  +  &tm2(m,u),  where 


brMi(u,  v)  —  /  (Vu  •  Vu  +  uvjdxdy  —  /  T^puvdi 

bTM2i.U,  V 


>Br 


=  (  —  l)  /  uvdxdy  —  /  (Tr  —  TR,p)uvd£ 

J Qr  J Br 

+  —dt'ya2  /  uvdi  +  a2  eruvdxdy. 

V  Js  Jcir 


(4.2.10) 

(4.2.11) 


Now  the  problem  becomes  finding  u  E  V  such  that 


brMi(u,v)  +  bTM2(u,v)  =  F(v)  Vu  G  V.  (4.2.12) 

The  motivation  for  splitting  up  the  sesquilinear  form  bpM(u,v )  (as  well  as  splitting 
up  the  Steklov-Poincare  operator,  Tr)  is  that  the  entire  term  itself  is  not  strictly 
coercive,  but  a  portion  of  it  (i.e.  6rmi  (t  v  j)  can  be  shown  to  be.  We  will  also  need 
to  show  that  bpM2(u,v)  is  a  compact  operator.  This  is  similar  to  the  idea  discussed 
by  Kress  in  [64]  of  decomposing  an  integral  operator  into  a  compact  operator  and  a 
operator  with  a  bounded  inverse.  He  notes  that  this  idea  goes  as  far  back  as  1919, 
and  he  notes  its  use  in  the  context  of  nonsmooth  boundaries  with  edges  and  corners. 
However,  we  will  also  see  that  this  decomposition  allows  us  to  show  that  bpMiu,  v) 
satisfies  a  Garding’s  inequality.  This  is  a  class  of  problems  where  the  sesquilinear 
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term  can  be  written  as  the  sum  of  a  coercive  term  and  a  compact  perturbation  [25]. 

At  this  point,  we  present  the  well-known  Lax-Milgram  Lemma  from  Cakoni  and 
Colton  in  [10],  which  will  be  used  to  prove  the  well-posedness  of  the  established 
problem: 

Theorem  4.2.1  (Lax-Milgram  Lemma)  Assume  that  b  is  a  bilinear  functional  map¬ 
ping  V  x  V  — *  C  for  which  there  exist  constants  Di,  D2  >  0  such  that 

\b(u,v)\  <  Di  1 1 M 1 1 y  ||u||y  \/u,v  G  V  (4.2.13) 

and 

b(u,u )  >  D-2  \\u\\2v  Vw  G  V.  (4.2.14) 

Then,  for  every  bounded  conjugate  linear  functional  F  :  V  — >  C  there  exists  a  unique 
element  liGf  such  that 

b(u,v)  —  F(v)  Vu  G  V.  (4.2.15) 

Furthermore,  ||«||y  <  D3  ||F||y/;  for  D3  >  0,  a  constant  independent  of  F. 

We  will  also  require  use  of  the  following  theorem  established  in  Cakoni  and  Colton 
in  [10]: 

Theorem  4.2.2  (Trace  Theorem)  Let  D  C  M2  be  a  simply  connected  bounded  domain 
with  dD  in  class  C2.  Then  there  exists  a  positive  constant  _D4  such  that 

IMIiP/^SD)  —  ^4  IMIffi(D)  ^ U  ^  H  (^)- 

Not  only  does  the  trace  theorem  state  the  trace  operator  can  be  “extended  as  a  contin¬ 
uous  mapping  from  F[l(D)  — »  Hll2(dDf%  but  also  “this  extension  has  a  continuous 
right  inverse”  [10].  That  means  that  for  any  /  G  FfL^2(dD),  there  exists  a  «G  H\D) 
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such  that  the  trace  of  u  is  /,  and  H^ll i^i(z?)  —  -^5  ll/ll.ffi/2(aDp  where  D5  >  0  is  a 
constant  independent  of  /  [10].  This  will  be  extremely  important  when  proving  the 
coercivity  and  continuity  properties  of  the  sesquilinear  term. 

We  are  now  ready  to  prove  our  main  theorem. 

Theorem  4.2.3  The  variational  problem  (4.2.1)  is  well-posed:  a  solution  u  G  V 
exists,  is  unique,  and  there  exists  a  C  >  0,  such  that 


lull  <  C 


+  us  +  u  +  u  +  sru 


Proof: 


We  first  must  show  that  &tmi(u,  v)  is  strictly  coercive,  i.e.  (4.2.14): 


bTMi(u,u)  —  /  (Vu  •  Vu  +  uujdxdy  —  /  T^puudi 

JO.R  JB  R 

1 2  ,  n  II  . .  II 2 


—  +  -^6  llMlljjl/2(Bfl) 

—  Il'ull_ff1(ni{)  +  -^7  IMI^Hr) 
>  (1  +  D7)  ||u||^i(r2i?)  • 
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We  also  show  that  1)tm\ (m,  v)  is  bounded,  or  continuous,  i.e.  (4.2.13): 


\bTMi(u,v)\  = 


(Vu  ■  Vu  +  uv\dxdy  —  /  Tptpuvdi 


< 


'  ^R 


V«  •  Vu  +  uv )  dxdy 


+ 


Tr  puvdd 


'Br 


/  uv  dxdy + 

/  TR'PUvdl 

1  qr 

Jbr 

—  II^MII  L2  Ill’ll  L2  +  IMI  L2  IMI L 2  T  ^8  IMIffi/2(Bfl)  ll^’ll/f1/2(Bi{) 

-  IMItf1^)  II 

=  (1  +  Dg)  ||m||^1(qh)  IMI/Z^Or)  • 


Thus,  as  in  Cakoni  and  Colton  in  [10],  we  may  apply  the  Lax-Milgram  Lemma  and 
state  there  exists  a  bijective  bounded  linear  operator  Bi  :  V  — >  V  with  bounded 
inverse  such  that 

bTMi(u,v)  =  {Bxu,v)  Vu  G  V. 

For  bpAniu,  v),  we  split  up  the  individual  terms  and  show  that  each  is  a  compact 
operator,  and  deduce  that  a  linear  combination  of  these  operators  is  compact  as  well. 
We  introduce  the  bounded  linear  operator  B2  :  V  — »  V: 


brM2(u,v)  =  ( B2u,v ) 

=  (  -  l)  (u,  v)L2{Qr)  -  {(Tr-  Tr,p)u ,  V )  +  (^5t~/a2)  (u,  v)L2{Ur) 

+  a2(£ru,v)L 2(f2fl)  Wv  e  V. 

Cakoni  and  Colton  in  [10]  prove  that  the  portion  of  the  operator  representing  (u,  v) l2(£ir) 
is  compact,  based  on  a  lemma  stating  the  injection  or  imbedding  of  H1{Qr)  into 
L2{Qr)  is  compact.  It  was  established  earlier  in  Theorem  3.4.8  that  Tr  —  Tptp  is  a 
compact  operator. 
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Following  Monk  in  [75],  we  define  the  space  L2  (Qr)  with  inner  product  as 


(■ u,v)L2r{nR)  =  (eru,v)  Vu,u  e  L2£r(nR). 


Then  the  norm  on  L2  (Q,r)  is  equivalent  to  the  standard  L2{VtR)  norm.  This  estab¬ 
lishes  that  the  term  a2(eru,  u)L2(nH)  is  compact.  Therefore,  we  can  deduce  that  the 
operator  B2,  a  linear  combination  of  compact  operators,  is  compact.  We  establish 
that  F(v)  is  bounded: 


\F(V)\  = 
< 


Jvdt  —  /  T Rusvdi  +  —dt'ya"  /  uvdi - /  uvdi  +  a  eruvdxdy 


>Br 


IBr 


V 


V  Js 


Jvdi 


IBr 


+ 

/  T Rusvd i 

+ 

—dt'ya2  /  uvdi 

+ 

—  [  uvdi 

Jbr 

V  Js 

V  Js 

a 2  /  eruvdxdy 


< 


< 


< 


f  ||  +  ||^i?Ms||  Ibll  +  C-z  ||u||  ||u||  +  CA  ||u||  ||u||  +  C5  ||£rw||  ||u|| 

dv}  ■  ~ 

~  Tru 1  ||u||  +  ||'hfi'Us||  ||u||  +  C3  ||u||  Hull  +  c4  IHI  ||u||  +  C5  ||eru||  \\v\ 

ch  /Ml  ~ 

"v\\  +  ||7rw*||  Hull  +  C2  ||tts||  \\v\\  +  c3  INI  INI  +  c4  ||u||  ||u|, 


dn 

+  C5  \\£ru\\  ||u|| 

<  Cl  ||«*||  ||u||  +  C2  ||uS||  IMI  +  C3  |N|  ||u||  +  C4  ||u||  ||u||  +  C5  ||£ru||  ||u 
=  Ci  I  ltd  1 1  +  C2  ||us||  +  C3  ||tt||  T  CA  I  lull  +  C5  H^rfi 


If  we  set  C  =  max{Ci,  C2,  C3,  C4,  C5},  then  we  have 


|F(u)|  <  C  ||u,:||  +  C  ||«-||  +  C  ||u||  +  C  p||  +  C  \\£ru\ 


=  C 


W  +  us  +  U  +  U  +  £rU 
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Since 


|F(u)|  —  &  11^*11  +  ll^s||  +  IMI  +  ||h||  +  ||£ru||  ||u|| , 

||F||  <  C  ||w*||  +  II^H  +  ||w||  +  ||u||  +  ||£r£t||  .  We  also  can  let  w  G  V  be  the  unique 
element  such  that 

F(v)  =  (w,v) 

which  is  true  due  to  the  Riesz  representation  theorem. 

At  this  point,  we  can  recast  (4.2.12)  as  hireling  u  €  V  such  that 

Biu  +  B2u  =  ( Bi  +  B2)u  =  w.  (4.2.16) 


Here,  we  cite  McLean  in  [73],  Theorem  2.33,  which  states  in  the  context  of  our 
problem  that  if  B  =  B\  +  B2,  where  B i  has  a  bounded  inverse  and  B2  is  compact, 
then  B  is  “Fredholm  with  zero  index,  and  hence  the  Fredholm  alternative  holds”  for 
the  equation  Bu  =  w.  As  a  result,  we  may  state  uniqueness  implies  existence. 

Thus,  in  proving  uniqueness,  brM(u,u)  =  0  implies 


We  further  assume  for  the  cavity  filling  material  9fJ(er)  >  0  and  A(er)  <  0,  whereas 
for  the  surface  S,  we  assume  fi  —  1  and  a  real-valued  rj  where  A(r/)  =  0  on  S. 

This  implies  that  the  imaginary  part  of  the  right-hand  side  of  (4.2.17)  is  less  than 
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or  equal  to  zero,  which,  in  turn,  implies 


9:(  /  TruuM )  <  0. 

Jbr 

Here,  we  apply  Theorem  3.2  presented  in  Ihlenburg  in  [54],  The  idea  of  this  theorem 

is  also  presented  in  several  different  sources  (see,  for  example,  Theorem  3.6  in  Cakoni 

and  Colton  in  [10],  Theorem  5.5  from  Angell  and  Kirsch  in  [6],  or  Theorem  2.12  from 

Colton  and  Kress  in  [19]).  This  theorem  states  that  a  solution  of  the  variational 

formulation  is  unique  if  $s(Gu,u)  <  0  holds  for  all  u  G  //1,/2(ra),  u  ^  0,  where 
du 

Gu  =  —  on  Ta,  where  Ta  is  an  artificial  boundary  enclosing  the  desired  domain. 
on 

Thus,  since 

A(  /  TriGuM)  <  0, 

Jbr 

then  u  =  0. 

The  application,  though,  is  in  the  context  of  the  exterior  boundary  value  problem 
(i.e.,  u  =  0  outside  of  f 1r).  Because  we  are  dealing  with  an  elliptic  operator  in  the 
form  of  the  modified  Helmholtz  operator,  the  solution  is  analytic.  Ihlenburg  further 
mentions  that  the  analytic  continuation  principle  (or  unique  continuation  principle  in 
Cakoni  and  Colton)  can  be  applied  to  state  u  =  0  in  the  interior  (i.e.  Dr).  Therefore, 
&tm(m,u)  =  0  implies  u  =  0,  so  the  solution  is  unique,  which  by  the  Fredholm 
alternative  implies  existence. 

Finally,  the  Fredholm  alternative  also  implies  the  boundedness  of  the  inverse  of 
(Hi  +  B-2),  which  for  constants  D0,  C  >  0  yields 


Ml  <  D0  ||F||  <  Do 

=  C 


c 


M  +  us  +  u  +  u  +  eru 


W  +  us  +  u  +  u  +  sru 
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This  completes  the  proof.  □ 

Thus,  we  have  established  well-posedness  of  the  solution:  existence,  uniqueness, 
and  an  estimate  which  bounds  the  solution.  From  this  estimate,  a  continuous  depen¬ 
dence  on  the  initial  (or  forcing)  data  can  be  determined. 


V.  Finite  Element  Analysis 


5.1  Construction  of  the  Discrete  Problem 

Now  we  transition  to  constructing  the  discrete  problem  from  the  variational  for¬ 
mulation  in  (4.2.1),  which  will  provide  the  foundation  to  solve  numerically.  Consider 
the  following  statement  from  Rjasanow  and  Steinbach  in  their  book  The  Fast  Solution 
of  Boundary  Integral  Equations : 


When  boundary  integral  equations  are  approximated  and  solved  numeri¬ 
cally  ,  the  study  of  stability  and  convergence  is  the  most  important  issue. 

The  most  popular  methods  are  the  Galerkin  methods  which  perfectly  fit 
to  the  variational  formulation  of  the  boundary  integral  equations.  [89] 

Not  only  does  Ihlenburg  in  his  book  Finite  Element  Analysis  of  Acoustic  Scattering 
add  this  is  an  area  of  ongoing  research,  but  he  also  describes  exactly  the  coupling 
approach  we  intend  to  use: 


In  practice,  one  is  not  necessarily  interested  in  computing  the  far-held 
results  directly  from  the  discrete  model.  Rather,  one  may  use  the  coupled 
finite-infinite  element  discretization  to  obtain  an  approximate  solution  of 
the  near  held  problem  in  (the  bounded  domain).  In  the  second  step,  one 
then  computes  the  far-held  pattern  from  the  Helmholtz  integral  equation, 
using  the  numerical  solution  on  a  “collection  surface”  in  the  near  held. 

If  such  an  approach  is  taken,  the  discretization  with  infinite  elements  is 
effectively  used  for  “mapping”  numerically  the  far-held  behavior  onto  the 
near  held.  [54] 

Applying  this  to  our  problem,  we  previously  derived  an  expression  for  the  helds 
exterior  to  the  cavity  in  (3.2.10),  and  also  obtained  a  variational  formulation  for  the 
near  helds  in  the  cavity  in  the  bounded  domain  (see  (4.2.1),  (4.2.8),  and  (4.2.9)). 
Once  the  solution  is  obtained  from  the  bounded  problem,  it  is  applied  to  obtain  the 
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desired  expression  of  the  far-held.  Clearly,  the  “collection  surface”  described  above 
corresponds  to  our  semicircular  artificial  boundary,  Br. 

The  finite  element  method  is  one  of  the  “most  powerful  methods  for  approximat¬ 
ing  solutions  to  PDEs”  [35].  It  is  based  on  the  fundamental  concept  of  writing  the 
boundary  value  problem  in  weak  or  variational  form,  as  in  (4.2.1),  and  using  the 
Galerkin  method  to  solve  the  equation  on  a  finite  dimensional  subspace  14  of  the 
solution  space  V,  resulting  in  a  linear,  finite  system  of  equations.  A  basis  of  piece- 
wise  polynomials  is  chosen  for  the  finite  dimensional  subspace  so  that  the  matrix  of 
the  linear  system  is  sparse.  We  will  consider  three  aspects  of  constructing  14:  the 
triangulation  of  the  domain,  working  with  functions  Vh  €  14  that  are  “piecewise  poly¬ 
nomials”,  and  creating  a  basis  in  the  space  14  whose  functions  have  small  supports 
that  are  easy  to  describe  [18]. 

For  the  first  aspect,  we  assume  the  bounded  domain  Hr  is  covered  by  a  family  of 
quasiuniform  triangular  subsets.  We  define  rh  =  {K}  as  the  partition  (or  triangula¬ 
tion  or  mesh)  of  Hr  where  each  K  e  t>(  represents  a  triangle  (i.e.  each  finite  element). 
These  finite  elements  form  an  exact  partition  of  Hr]  that  is,  Hr  =  U  KerhI<  [7]. 

For  an  arbitrary  triangle  K,  we  denote 

hK  =  diam(A')  =  max{||x  —  y||  j  x,  y  G  K}, 

which  corresponds  to  the  length  of  the  longest  side.  We  further  define  the  mesh  size, 
h,  of  the  partition  r^,  as 

h  =  max  h  k  , 

K&rh 

and  Sk  as  the  diameter  of  the  largest  circle  inscribed  in  K  [7].  A  family  of  triangula¬ 
tions  is  regular  if  there  exists  a  constant  c0  such  that  hKj  Sr  <  c0  for  all  K  and  h  — >  0 
[7].  We  further  note  from  Atkinson  and  Han  in  [7]  that  a  family  of  triangulations 
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{tTj}  is  quasiuniform  if  each  triangulation  Ty  is  regular,  and  there  exists  a  constant  c0 
such  that 

min  Hk/  max  hj<  >  Co- 

For  the  second  and  third  aspects  in  constructing  14,  following  [104],  we  have  the 
finite-dimensional  subspace  14  of  the  test  space  V : 

Vh,  =  {vh  £  H\nR)  :  Vh | k  £  Pi,K  G  r^}, 
where  {4>*  (a?) is  a  linear  nodal  basis  of  14-  Each  Vh  G  14  is  expressed  as 

N 

I’h  = 

3=1 

(That  is,  each  i>h  can  be  written  as  a  linear  combination  of  the  basis  vectors.)  Ac¬ 
cording  to  Ciarlet  in  [18] ,  Vh  G  14  is  the  key  to  all  convergence  results,  and  helps  yield 
simple  computations  of  coefficients.  We  also  note  14  is  closed  in  V  and  14  — >•  V  as 
h->  0. 

Therefore,  having  constructed  14,  the  fully  discrete  problem  is  to  find  ruvh  G  14 ,  n  — 
1,2,...,  iV,  such  that 


bTM«,vh)  =  Fn(vh)  Vvh  G  14,  (5.1.1) 

where  i’h)  and  Fn(yh)  are  dehned  as  in  (4.2.8)  and  (4.2.9),  respectively.  The 

bilinear  form  6tm('«4  vh)  can  be  written  as  the  sum  of  a  positive  definite  (coercive)  op¬ 
erator  and  a  compact  operator,  as  seen  previously  for  (4.2.12).  According  to  Pomp  in 
[83]  and  Hildebrandt  and  Wienholtz  in  [44] ,  this  is  equivalent  to  satisfying  a  Garding’s 
inequality;  thus,  the  discrete  problem  (5.1.1)  has  a  unique  solution  and  Cea’s  lemma 
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applies: 


\\un  —  v%\\v  <  C  inf  \\un  —  Vh\\v  ■  (5.1.2) 

Vh&Vh 

Furthermore,  “for  a  mesh  width  h0  >  0  and  a  constant  C  >  0  such  that  for  all  h  <  ho, 
the  condition  of  uniform  14  ellipticity  applies” : 

\bTM{vh,vh)\  >  c  \\vh\\2  \/vh  e  Vh, 

which  means  the  linear  system  has  a  unique  solution  due  to  the  Lax-Milgram  lemma, 
for  all  h  <  h0  [44], 

This  is  considered  the  /i- version  of  the  hnite  element  method  (and  why  we  write 
Vh),  where  we  achieve  convergence  by  refining  the  mesh  (h  — y  0).  Another  option 
is  to  increase  the  polynomial  degree,  referred  to  as  the  />- version  of  FEM.  It  is  also 
possible  to  combine  these  methods  through  an  hp- version.  Nevertheless,  the  choice 
of  method  depends  on  the  knowledge  of  regularity  of  the  problem  [7].  In  addition, 
since  Vh  C  V  for  each  h,  this  approximation  is  conforming  [75]. 

5.2  Error  Analysis 

Recall  Cea’s  lemma,  (5.1.2),  which  simply  shows  that  to  estimate  the  hnite  solution 
error  (or  discretization  error)  \\un  —  r4||y,  we  only  need  to  estimate  the  approximation 
error  ||un  —  Vh\\ v-  The  inhmum  portion  of  Cea’s  lemma  “characterizes  the  exact 
solution  in  the  subspace  that  is  spanned  by  the  FE  shape  functions.  If  the  inhmum 
is  reached  on  an  element  i’h  G  14,  then  this  element  is  the  best  approximation  of 
un  in  the  subspace  14”  [54].  We  further  note  that  Cea’s  lemma  is  a  quasioptimal 
error  estimate  because  up  to  constant  C,  the  discretization  error  is  bounded  by  the 
approximation  error.  An  optimal  estimate  would  be  C  =  1  [75]. 

It  is  clear  the  Galerkin  method  is  tied  directly  to  the  H 1  or  energy  norm,  so  we 
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would  expect  the  best  approximation  to  the  solution  to  be  measured  in  this  norm 
[36].  H1  error  norm  is  “relevant  for  exterior  problems  if  one  is  interested  in  the 
far-held  response”  [54],  Recall  that  we  are  computing  the  scattered  held  from  a 
boundary  integral  equation  in  the  exterior  of  the  cavity,  where  we  will  use  the  finite 
element  data  on  the  artificial  boundary  (or  collecting  surface),  Br.  Ihlenburg  adds 
that  “since  this  integral  equation  involves  u  and  its  normal  derivative,  the  H 1  error 
norm  is  appropriate  for  error  control  in  the  near  held.  The  situation  is  different  for 
interior  problems.  Here,  the  L2  error  norm  may  be  more  appropriate  for  error  control” 
[54]. 

However,  the  L2  space  is  not  always  well-suited  for  approximating  functions  in 
regions  which  are  rapidly  oscillating  with  small  amplitudes.  Small  details  may  be  lost, 
so  H1  incorporates  not  only  the  difference  in  the  functions  but  also  in  the  gradients 
[38].  Therefore,  in  the  context  of  the  solution  to  the  interior  problem,  our  goal  is  to 
quantify  the  finite  solution  error  ||un  —  ii)[  1 1  in  both  the  L 2  and  H 1  norms. 

In  general,  there  are  two  types  of  error  estimation:  a  priori  and  a  posteriori. 
An  a  priori  estimate  is  the  “error  to  be  expected  in  a  computation  to  be  done,” 
whereas  a  posteriori  estimates  are  “generated  in  the  course  of  computation”  [87]. 
More  specifically,  according  to  Ihlenburg,  a  priori  error  estimation  is  descibed  as 
follows  : 


The  error  function  is  estimated  in  a  suitable  functional  norm  without  quan¬ 
titative  input  from  the  computed  solution.  The  estimates  are  based  on  the 
approximation  properties  of  the  subspace  where  the  numerical  solution  is 
sought  and  on  the  stability  properties  of  the  differential  operator  or  vari¬ 
ational  form.  The  estimates  are  generally  global;  i.e.,  the  error  function  is 
estimated  in  an  integral  norm  computed  over  the  whole  solution  domain. 
[54] 


He  adds  that  a  posteriori  error  estimation  can  be  described  as  follows: 
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The  error  function  is  estimated  employing  the  computed  solution  of  the 
discrete  model  as  data  for  the  estimates.  In  practice,  the  estimates  are 
usually  part  of  an  adaptive  mesh  refinement  methodology.  For  mesh  re¬ 
finement,  one  needs  local  information  on  the  error.  A  posteriori  error 
estimation  should  therefore  be  given  in  a  norm  that  is  defined  on  a  single 
element  or  on  a  patch  of  adjacent  elements.  [54] 


Furthermore,  the  a  priori  error  bounds  are  asymptotic,  not  absolute.  That  is, 
“the  bounds  will  not  tell  how  small  the  error  is  when  the  solution  is  approximated  on 
a  particular  mesh.  Instead,  the  bounds  show  how  the  error  decreases  as  the  mesh  is 
refined”  [36].  Accordingly,  local  mesh  adaptation  is  not  used.  This  is  in  contrast  to 
a  posteriori  error  bounds,  which  are  expressed  in  terms  of  residuals  of  approximate 
solution,  not  in  terms  of  powers  of  a  mesh  size  or  constants  depending  on  the  exact 
solution.  [87].  Therefore,  we  are  providing  only  a  priori  error  estimates. 

Up  to  now,  we  have  not  determined  the  symmetry  of  the  sesquilinear  term  bTM(u,  v ). 
It  has  not  been  a  factor  for  the  existence  and  uniqueness  proofs,  so  in  order  to  estab¬ 
lish  the  error  bounds  for  the  L 2  and  H 1  norms,  we  will  follow  [92]  and  assume  the 
more  general  case  that  v)  is  nonsymmetric,  but  is  still  bounded  as  before. 

In  addition,  the  following  theorem,  which  is  similar  to  one  presented  by  Van  and 
Wood  in  [105],  will  determine  uniform  convergence  estimates: 

Theorem  5.2.1  Let  un  G  V  and  vf  G  14  be  the  solutions  to  (4.2.1)  and  (5.1.1), 
respectively,  for  Fn  G  V' .  Given  e  >  0,  there  exists  an  h0  =  h0(e)  such  that  for  all 
0  <  h  <  ho,  then 

llM'  —  uh  llz,2(nfl)  —  6  \\u  ~  uh\\v  ■  (5.2.1) 

Furthermore,  if  given  e  >  0,  there  exists  an  hi  =  hi(e)  such  that  for  all  0  <  h  <  hi, 
then 

K  <  Ce  ||F"||I1(njt)  (5,2,2) 
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for  some  positive  constant  C  independent  of  h.  It  follows  that 


ll“”-<llLW<C'62||f’”llil(n„).  (5.2.3) 

Also  we  will  use  the  following  Lemmas,  the  first  of  which  is  proven  in  both  [92] 
and  [105],  and  the  second  is  proven  in  [92]: 

Lemma  5.2.2  Let  D  —  {/  :  /  G  L2( flR),  ||/||L2(Qfl)  =  1}  be  the  unit  sphere  in 
L2( LIr).  For  f  G  D,  let  W  be  the  set  of  solutions  w  G  V  to  the  auxiliary  /  adjoint 
bilinear  form 

bTM{w,v)  =  (f,v)  Vu  G  V 

Then  W  is  compact  in  V . 

Lemma  5.2.3  Let  V  be  a  fixed  compact  subset  of  Then  given  any  e  >  0, 

there  exists  an  h0  =  h0(e,V)  such  that  for  each  u  G  V  and  each  0  <  h  <  ho,  there 
exists  a  Vh  G  14  satisfying 

|| u  -  vh\\v  <  e. 

Proof  of  Theorem  5.2.1: 

We  use  an  Anbin-Nitsche  duality  argument,  assuming  brM(u,v )  is  nonsymmetric, 
and  follow  the  reasoning  in  [92]  and  [105].  Let  ')  be  the  adjoint  bilinear  form 

to  •)  defined  by 


b*TM(u,  v )  =  bT M(v,  u)  Vu,  v  G  V. 


We  know 

b*TM(w*,v)  =  (g,v)  Vu  G  V 
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has  a  unique  solution  for  all  g  e  L2(flR),  and 


\w  llu  —  c \\y\\L2(nR)  ■ 


If  we  view  un  -  as  a  linear  functional  in  L2(QR),  then 


I  un  -  un 


h\\L2(QR)  =  SUP  («”  -<,£>■ 
Hfflli2(nH)=l 


Then,  for  Vh  G  14  ,  we  have 

(g,un-unh)  =  (un-unh,g) 

=  b*TM{w*  ,un  —  v%) 

=  b*TM(w*  —  Vh,  un  —  Ufr)  +  b*TM(yh,  un  —  u%) 

=  bTM{u11  —  11%,  w*  —  i>h)  +  bTM(un  —  u1^,  i>h) 

=  bTM(un  ~  <,  w*  -  vh) 

<  c\\un  -<||y  |K  ~vh\\v. 

By  using  an  approximation  property  argument  from  Strang  and  Fix  in  [99], 
can  choose  Vh  such  that  ||w*  —  Vh\\v  <  e  ||tc*||y.  Therefore, 


(un  ~  K,  g)  <c  || un  -  <||y  e  \\w* 


Iv 


<Cl  ||«  -<||y6C2||^||L2(ni?) 

in  n  \  ^  1 1  77  to  1 1 

sup  {u  — uh,g)<  sup  ci  ||m  —  uh\\vec2 


Ils,llz,2(r2^j)  — 


=1 


K  -  <llL2(n*)  ^  Ce^u 


n  un 


h\\V  ' 


L2^R) 


we 
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which  shows  (5.2.1).  To  show  the  second  estimate  (5.2.2),  as  in  [92]  and  [105],  we  set 


pn  _ 


un  = 


u 


ul  = 


u 


h 


l£2(fij0 


lL2(Ofl) 


Then  we  have  the  corresponding  variational  problems 


bTM(un,v)  =  Fn(v)  WeV, 

bTM(uh,Vh)  =  Fn(vh)  Vi’h  E  14, 


for  which  we  apply  Cea’s  lemma, 


\un-ul\ 


v  <  C  inf 

v  vhevh 


| un  -  vh\ 


v  ■ 


From  Lemma  5.2.2,  the  set  of  solutions  for  &TM(hn,  v)  =  Fn{y ), 


L2(Oh) 


is  compact  in  V,  thus  we  can  apply  the  density  argument  in  Lemma  5.2.3  to  get 


1, 


inf  | \un-vh\\v  <  e 

vh& 4 


for  0  <  h  <  h0(e,un).  Thus  we  conclude  that 


\\un  -  Unh\\v  <  Ce. 

1 1  7/n  _  1 1 

Since  \\un  —  u^Wv  —  t  ^  n — then 

II  n  1 1 V  Pn 

IT  \\L^(Qr) 

\\un-ul\\v<Ce\\Fn\\L2{QR)n 


In  summary,  we  are  generating  a  priori  error  estimates,  and  are  simply  applying 
previously  established  theory  to  confirm  the  error  bounds  for  our  problem. 


77 


5.3  Stability  Analysis 


For  stability  analysis,  we  initially  follow  Van  and  Wood  in  [105]  precisely  and  ex¬ 
press  the  Newmark  scheme  in  a  three-step  formulation.  We  start  with  the  discretized 
PDE: 

—A  un+2  +  a2£run+ 2  =  a2£run+ 2 

=  a2£r  un+1  +  5tiin+l  +  ®  (1  -  2 j3)iln+l  . 

Using  (3.1.7)  to  remove  un+1  and  (3.1.4)  to  remove  ■un+1,  we  obtain 

-A  un+2  +  a2£run+ 2 

=  a2£r  [>+1  +  5tun  +  (5t)2(  1  -  7 )un  +  (At)2 7nn+1  +  ^(1  -  2 :/3)un+1  . 

Using  (3.1.3)  to  remove  un ,  followed  by  using  (3.1.6)  to  remove  iln+l  and  un+1  sepa¬ 
rately,  and  finally  applying  (3.1.5)  to  remove  un+1  and  un  separately,  we  obtain 

(3Aun+2  +  (^  -  2 0  +  7)A un+1  +  (^  +  p  -  7)A un  =  /3a2£r(un+ 2  -  2 un+1  +  un ). 

Therefore,  the  final  variational  form  of  this  equation,  fully  discretized  in  both  space 
and  time,  using  v%  G  14,  is 

^4«+2-2<+1  +<),»») 

+  a(/3Au”+2  +  (i  -  2/3  +  7)  A<+1  +  (!  +  /?-  7)  A<,  vh)  (5-3.1) 

=  /3Gn+2(vk)  +  (1  -  2/3  +  'y)Gn+1(vh)  +  (\+P~  l)Gn(vh)  V17,  £  14- 
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This  is  an  expected  result  in  the  literature,  as  seen  in  Van  and  Wood  in  [105],  Quar- 
teroni  and  Valli  in  [85],  and  Allaire  in  [4],  We  further  define  for  our  problem 

a(it%,Vh)  —  /  Viih  ■ 'Vvhdxdy  —  /  TRvU^Vhdt  +  —St^a2  /  u^Vhdt 
Jclr  Jbr  V  Js 

and 


Gn(u)  =  /  Jnvhdi  —  /  VRu8’nvhd£  +  -5t^a2  /  unvhd£  -  -  /  ~unvd£. 
Jbr  Jbr  V  Js  V  Js 


We  note  that  a(u%,  Vh)  is  symmetric,  but  not  strictly  coercive.  Therefore,  as  in  the 
previous  section,  we  rewrite  a(u%,  Vh)  =  a i(v%,  Vh)  +  a 2(u^,  Vh)  as  the  sum  of  a  coercive 
term, 

ai(ul,vh)  =  /  V <  ■  Vvhdxdy  -  /  TR,Pulvhd£ 

JflR  JBR 

and  a  compact  term, 


a2«Wh)  =  ~  [  (' Tr  -  TR,P)ulvhd£  +  -dt'ya2  [  u\vhd£. 

Jbr  V  Js 

Therefore,  as  cited  earlier  in  Pomp  in  [83],  as  well  as  Lemma  22.38  in  Zeidler  in  [117], 
a(u%,  Vh)  satisfies  a  Garding’s  inequality.  That  is,  there  exists  two  positive  constants 
k  >  0  (Zeidler  says  k6R)  and  v  >  0  such  that 

a(u,u)  +  k,\\u\\2L2  >  v  \\u\\y  WuEV.  (5.3.2) 

Therefore,  the  term  a(v%,  Vh)  is  symmetric  and  satisfies  a  Garding’s  inequality. 
Therefore,  according  to  Allaire  in  [4],  the  eigenvalue  problem 


a(wh,  vh)  =  A h{wh,  vh)  Mvh  e  Vh 
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has  a  finite  number  of  eigenvalues  satisfying 


—  K  <  Aft,!  <  Xh,2  <  •  •  •  <  A h,M  <  oo, 

where  dim  14  =  M  for  corresponding  orthonormal  eigenvectors  11^4,11^,2;*  '•  •  •  ,u>h ,m- 
As  in  [105],  writing  =  w^i  and  A;  =  A/^,  and  replacing  vj,  for  ry  in  (5.3.1),  we 
have 

a(u%,Wi)  =  a(wi,  u'l)  =  \{wi,v%)  =  \{u^,Wi). 

This  is  also  equivalent  to,  in  referring  to  (5.3.2),  as 


a(v%,Wi)  +  K(u%,Wi)  =  \i{u%,Wi)  +  K{v%,Wi)  =  (Aj  +  K){v%,Wi) 


and 

^(Sr«+2-2<+‘  +  u;),Wi) 

+  (A.  +  K)  (/3A«”+2  +  (1  -  2/?  +  7)A<+1  +  (1  +  /?  -  7)  A<, «,.) 

=  PGn+2{Wi)  +  (l--2(3  +  7 )Gn+\wi)  +  +  (3  -  7 )GnK) 

=  G*. 

Without  loss  of  generality,  we  further  assume  er  —  1  and  G*  =  0,  and  from  [105],  we 
have  the  eigenvalue  equation:  find  Xh  and  Uh  G  14  such  that 

a(uh,  vh)  +  n(uh,vh)  =  (Ah  +  K)(uh,vh)L2r(nR)  Vty  G  14.  (5.3.4) 

Now  we  have  that  a(-,-)  +  «;(•,•)  is  coercive,  so  we  have  Aj  +  k  >  0,  where  the 
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corresponding  orthonormal  eigenvectors  satisfy 


(wi  >  Wj  )  Llr  (f2fl)  &ij 


Hence,  by  replacing  u £  =  YmLi  u7wh,i  into  (5.3.3),  we  get  for  i  =  1,  2, ...,  M, 


1  f?,n+2  _  2?/n+1  4-  7/"^ 

(«)d  •  •  +  ') 

+  (A.  +  k)(/3A<+2  +  (i  -  2/3  +  7)A<+1  +  (i  +/3  -  7>A<)  =  0. 


(5.3.5) 


We  rewrite  this  three-level  scheme  as 


_  2  —  (Aj  +  K)(5t)2(\  —  2(3  +  7)  n+1  1  +  (Aj  +  K)(6t)2(\  +  /3  —  7)  n 
"  l  +  (A,  +  «)(5f)2^  l  +  (A.8  +  «)(<5t)2^  Ui 

=  bnUi+1  - 


Thus,  (5.3.5)  can  be  written  in  matrix  form  as 


<+2 

611 

—bi2 

<+i 

=  3, 

<+' 

<+1 

1 

0 

< 

< 

where  Bi  denotes  the  iteration  matrix,  and  for  stability,  the  Von  Neumann  necessary 
condition  is  p(Bi)  <  1  where  the  spectral  radius  p(Bi)  denotes  the  maximum  modulus 
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of  the  eigenvalues  of  Bj.  Therefore,  we  consider  the  characteristic  equation  of  Bf 


{C  0  &n  — b\2 

— 

0  cj  [i  0 


—  (2  —  (bn  +  b\-2  —  0. 


0 


and  calculate  the  roots  of  this  equation  which  are  the  eigenvalues  of  the  matrix  B{. 
The  roots  are 

In  order  to  analyze  when  |Ci|  ,  IC2I  <  1,  we  must  calculate  the  discriminant,  A: 


A  —  b\x  —  4&i2 

2  —  (A  i  +  K,)(St)2(^  —  2/3  +  7)  1  +  (A*  +  n)(8t)2(^  +  /3  —  7) 

1  +  (Aj  +  K)(5t)2(5  1  +  (A  i  +  K,)(St)2/3 

-4(Aj  +  n)(8t)2  +  (Aj  +  K)2(8t)4[( 7  +  0.5)2  -  4/3] 

[1  +  (A  i  + n)(8t)2f3}2 

This  expression  is  confirmed  in  Allaire  in  [4].  Since  the  denominator  is  always  positive, 
in  analyzing  the  sign  of  A,  without  loss  of  generality,  we  will  assume  henceforth  that 

A  =  — 4(Aj  +  K)(5t)2  +  (Aj  +  (7  +  -)  —  4/3  . 

Consider,  then,  the  following  two  cases:  Case  1:  A  <  0 
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This  implies  that  the  roots  Ci>C2  are  complex  conjugates,  so  we  require 


161  =  161=  bf  + 


%  -  4&! 


bj  i  b2u  -  4  h 


Thus,  if  &i2  <  1,  then 


_  1  +  (A*  +  K)(5t)2(\  +  /3  —  7) 

12  “  l  +  (A,,  +  «)(^)2/3  -  ’ 


which  is  equivalent  to  ^  +  B  —  7  <  B,  or  7  >  As  a  result,  we  have  that  for  this 

z  2 


(A*  +  K)(5t)2  (q+^)2-4 /3  <4  =»  7 


Case  A-  A  >  0 

This  is 


(Aj  +  n){5t)2  (7  +  -)  —  4/3  >  4. 


Then,  without  loss  of  generality,  we  assume  Ci  <  (2,  so  we  require 


—  1  ^  Ci  —  —  — -  7  —  "T - —  C2  ^  1* 

-  Si  2  2  2  2 


For  —  1  <  Ci,  we  have 


611  v^A  ,  ,  ,  .  , 

—  1<— - —  =7  1  +  bi2  >  —bn. 


Whereas  for  C2  <  1  we  have 


611  +  v/A 


<  1  =>  1  +  bu  >  bn- 


These  together  imply  that  1  +  &12  >  |&n| .  However,  we  see  that  this  results  in  two 
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inequalities  that  must  be  satisfied  [105].  For  1  +  612  >  — &11,  we  have 

1  7 

(A t  +  n)(6ty  ~  2~  13 ’ 

but  for  1  +  bi2  >  &n  we  have 

(Ai  +  K)(5t)2  -  2  ;  4 

Therefore,  in  order  to  satisfy  both  inequalities,  we  can  reduce  this  to 

1  7 

- >  -  -  8. 

(A,  +  «)(5t)2  -  2  ; 

To  summarize,  for  Case  1  (A  <  0),  we  have 

(Xi  +  K)(St)2  (7  +  i) 2  —  4/3  <4  implies  7  >  i. 

Whereas  for  Case  2  (A  >  0)  we  have 

(A*  +  hi)(St)2  (7+i)2-4/3  >4  implies  (A*  +  K)(5t)2{2^  -  4/3)  <  4. 

We  further  assume  that  /3  >  0  and  7  >  -.  In  fact,  it  is  shown  in  Dautray  and 
Lions  that  the  scheme  is  unstable  for  7  <  -  [24],  We  first  observe  the  property 
that  (A;  +  hi)(6t)2  >  0.  Thus,  if  (7  +  |)2  —  4 8  <  0,  then  only  Case  1  applies. 
If  (7  +  7  —  4/3  >  0,  then  we  have  stability  as  long  as  (2 7  —  4/3)  <  0.  These 
statements  can  be  combined  to  deduce  that  we  have  unconditional  stability  for  7  <  2/3 
for  arbitrary  5t.  However,  if  7  >  2/3,  then  we  require 

4 

max(A;  +  hi) (St)2  <  (27_4/3), 
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which  is  a  conditional  situation  that  depends  on  finding  max,  A*  and  k. 

It  can  also  be  shown  that  when  7  =  -,  a  conditional  stability  arises  dependent  on 
finding  the  values  of  \  and  k.  Also,  Van  and  Wood  state  that  7  =  -  is  not  always 
the  best  value  to  use,  and  that  in  practice,  a  value  of  7  >  -  “is  often  used  to  damp 
out  the  higher  frequencies  while  preserving  the  more  accurate  lower  ones”  [105]. 
Therefore,  we  can  summarize  these  results  in  a  theorem: 

Theorem  5.3.1  The  Newmark  scheme  for  the  TM  variational  problem  is  uncondi¬ 
tionally  stable  for  arbitrary  St  >  0  satisfying 

2/3  >  7  > 

As  with  the  error  analysis  section,  during  our  stability  analysis  we  have  built  upon 
existing  theory  and  analysis,  and  refined  it  to  apply  to  our  problem. 
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VI.  Numerical  Simulation 


6.1  Finite  Element  Approximation 

The  goal  in  this  section  is  to  set  the  framework  for  the  numerical  implementa¬ 
tion.  Starting  with  the  triangulation  and  basis  as  described  earlier,  the  fully  discrete 
problem  is  to  find  v%+1  G  14,  n  =  1,  2, . . . ,  N,  such  that 

bTM(unh+1,vh)  =  Fn+1(vh)  Vvh  G  14,  (6.1.1) 


where  6™(^+V'/i)  and  Fn+1(vh)  are  defined  as  in  (4.2.8)  and  (4.2.9),  respectively. 
Since  each  Vh  G  14  and  u )(+1  G  14  is  expressed  as 


Vh 


and 


u't+1  = 


N 

£<+1< 

i— 1 


,h , 

l  5 


then 

N  N  N 

*— 1  1=1  1=1 

Due  to  linearity,  this  yields 

N  N  N 

£  «,!»«/(£  <+Vt  0?)  =  £  <^"+1(0h 

1=1  *=l  1=1 


In  order  for  equality  to  hold  for  all  Vj,  it  follows  that 


N 


W£<+V!‘,0-)  =  i?”+1(0* )  forj  =  1.....JV. 


Applying  linearity  again  yields 


N 

Y,  #K+1  =  Fn+1(<pb)  for  j  =  1, ....  N. 

1=1 

Then,  referring  to  (4.2.8)  to  discretize  brM(<f>i,<f>j)  ,  we  obtain  the  matrix  equation: 

N 

]T([A y  +  a2  [M\ji  +  ^StWlPh  +  [QhH+1  =  Fj+1  for  3  =  !,  N.  (6.1.2) 

i— 1  ^ 

Here  we  define  the  stiffness  matrix 


[K\ji=  /  Vtf-V&dxdy,  (6.1.3) 

JnR 

the  mass  matrix 

[M]ji  =  [  ertftfdxdy, 

JuR 

the  boundary  matrix 

\P\ji  =  f 
Js 

and  the  artificial  boundary  matrix 

[Qh  =  ~  [  Tntf^de.  (6.1.6) 

Jbr 

The  right-hand  side  vector  F™+1  is 


(6.1.4) 

(6.1.5) 


r  fhii>n+ 1 


^Rus’n+1^de 


—5t^a2 

V 


u 


n+lih 


d 


Br 

«+i 


3 


^  /  iin+  wdl  +  a 
V  Js 


3 

2  /  £run+1(j)^dxdy. 


(6.1.7) 


While  the  stiffness,  mass,  and  boundary  matrix  expressions  are  straightforward 
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for  the  left-hand  side  of  the  equation,  the  artificial  boundary  matrix  is  more  difficult 
to  evaluate.  Here,  we  refer  to  Steinbach’s  work  in  which  he  addresses  “defining  a 
computable  approximation  of  the  Steklov-Poincare  operator”  [96].  The  advantage 
of  his  approach  is  that  it  is  a  symmetric  approximation  that  does  not  involve  the 
hypersingular  operator.  Referring  back  to  (3.4.15),  we  have 

Tr  =  S~lZS~l 

where  Z  =  (^/  +  D)S  and  Z  :  H~1/2{Br)  — >  H1^2( Br ).  We  further  define  for  < p(r ')  = 
du(r') 
dnr, 

(Z(p)(r)  —  ~  [  G(r\r')(p(r')d9'  +  f  ^G(r\r  )  I  £(r,|s,Ws/)d0',d0/.  (6.1.8) 
2  Jbr  Jbr  dnr,  JBr 

Thus,  Z  admits  a  direct  Galerkin  discretization. 

We  also  need  to  further  define  Xh  and  Yh,  where  Xh  is  the  space  of  traces  of  finite 
elements  Vh  G  14  C  V  =  H1(BIr)  and  W  is  its  dual  space.  Thus,  we  have  the  spaces 
of  trial  functions: 


Xh  =  span{0fe}f=1  C  H1/2(Br)i 
Yh  =  span{^fe}^=1  C  H~1/2(Br). 

We  see  that  Xh  is  the  approximation  of  u\Br,  consisting  of  continuous  piecewise 
linear  functions,  and  Yh  is  the  approximation  of  dnu  \Br,  consisting  of  constant  basis 
functions.  Furthermore,  the  “piecewise  constant  boundary  functions  can  represent 
the  Neumann  derivatives  of  the  piecewise  linear  functions  exactly”  [45]. 


Therefore,  in  order  to  evaluate 


IQlji  =  = 


T T  c-1 
1h 


Zhs^ih, 


we  can  define 


[sh]i,k  =  (6.1.9) 

\zh\k  =  l^Br),  (6.1.10) 

[4]*,f  =  (6.1.11) 


Therefore,  the  artificial  boundary  matrix  and  consequently  the  left-hand  side  admits 
a  Galerkin  discretization. 

We  now  describe  how  to  assemble  and  compute  the  single  layer  potential  matrix, 
S.  The  technique  described  can  also  be  applied  in  constructing  Z .  We  follow  the 
technique  of  Van  and  Wood  in  [104]  using  a  standard  quadrature  with  the  midpoint 
formula: 


I&ks  =  <'Sv£d‘W»> 


=  Vi 


G(r\r')'ijj^(r')dr'dr. 


If  we  divide  the  segments  on  the  artificial  boundary  as  Bj.  =  (r*^,  t*^_|_1)  and  Bj  = 

(rfrri+i)i  we  then  further  define  ^  =  A'+1^ - -  and  ^  =  ,~+1^ - -.  Now,  applying 

the  quadrature  using  this  midpoint  formula,  and  noting  that  the  basis  functions  are 


piecewise  constant  over  the  segment  with  a  value  of  1,  for  Bk  ^  Bj  we  have 


G(r\r')^(r')dr'dr  «  G(f* |fj). 


Now  considering  the  case  along  the  diagonal  of  the  matrix  S,  or  when  Bk  =  Bf,  or 
=  we  wi^  have  singularities  for  the  first  two  terms  in  the  Green’s  function 
(3.3.7).  We  will  address  this  as  Van  and  Wood  did  in  [104]  and  have 


r')^- (r'Wr'dr  ~  /  K0(T)dr. 

a  I  n 


and  the  remaining  terms  of  the  Green’s  function  in  (3.3.7)  can  be  evaluated  as  de¬ 
scribed  in  Section  3.5.  Also,  the  above  expression  can  be  evaluated  exactly  as  provided 
in  [1]  for  the  term 

Kq{t  )dr . 

Therefore,  we  have  S'  as  a  square,  dense  boundary  element  matrix  of  size  /  x  k. 
This  matrix  is  invertible  and  each  entry  element  will  have  a  real  and  imaginary 
portion,  if  we  include  the  third  term  of  the  Green’s  function  in  (3.3.7).  The  same 
holds  for  the  matrices  Z  and  I.  These  three  matrices  will  ultimately  determine  the 
artificial  boundary  matrix,  Q,  as  described  in  (6.1).  Because  the  dimensions  of  Q 
are  not  the  same  as  the  stiffness,  mass,  and  boundary  matrices,  we  can  standardize 
it  to  the  other  dimensions  using  a  process  described  in  detail  by  Jin  in  [55].  The 
difference  between  the  artificial  boundary  matrix  and  the  boundary  matrix  is  that 
the  artificial  boundary  matrix  reflects  the  fact  that  each  boundary  element  interacts 
with  every  other  boundary  element,  hence  the  dense  square  matrix.  The  boundary 
matrix,  however,  reflects  the  fact  that  each  boundary  element  interacts  only  with 
adjacent  boundary  elements. 
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In  analyzing  the  right-hand  side  of  the  equation  given  by  (6.1.7),  we  first  note  as 


in  [104]  that  un+l  ~ 


so 


a  /  eru 

J 


N  „  N 

n+1(f)^dxdy  ~  -un+1  /  Erfy^dxdy  =  sy^\M]jiu^+l .  (6.1.12) 

i=i  i=i 


Similarly,  we  have  iin+l  pe  ,  which  implies  we  also  have  the  terms 


»  N  „  N 

—dt^a2  /  un+1<p^d^  ~  —dt^/a2  un+1  /  (f^cf^dl  —  —dt^a2's^[P]jiu7l+l.  (6.1.13) 
‘d  d  s  7  ■ _ i  d  s  7 
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«  «  E  *?+1  /  = -- 


n+1 


2=1 


Likewise,  we  can  discretize  the  following  term  as: 


(6.1.14) 


^inc,n+l 

dn 


)$d£  = 


BRj  |  (9umc’n+1(^-,  (n  +  l)hf) 
2  <9n 


(6.1.15) 


where  \Br.  |  is  the  length  of  the  discretized  segment  on  the  artificial  boundary, 
is  the  midpoint  of  the  discretized  segment,  and  (n  +  1  )5t  is  the  time  of  evaluation. 
These  also  apply  in  discretizing  the  following  term  as: 


TrU 
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inc,n+l  ^pli  d£ 
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Rj  1  ry  inc,n-\- 1 


2 

\B* 


B 


Rn 


U 


1 Br 


Th u*”’”+1(&,  (n  +  l)St) 

i{S-1(l/  +  D)«~’»+1K,',(n+l)*)} 

l{S-1[i«‘“’”+1K„(n  +  l)«) 

,«+Urrf_mndr,]} 


-irt 


2U 


inc,n+ 1 


(€j,  (n  +  l)St) 


E  I»hJ  (n  +  !)«)*]  }■ 


(6.1.16) 
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We  further  note  in  the  previous  cases  that 


dumc,n+ 1(£^  _|_  X)ht) 


inc,n-\- 1 


dn 


dGfa  |fe) 

dn' 


=  Vw 


=  VG-n&, 


for  outward  unit  normal  vectors  to  the  boundary  he  and  he, . 

Jbr 


Br\ 


IBr 


'FRhs’n+1(r)<^d£«  L^^Rfis’n+1(^), 


where 


=  *S'_1|tt2  J I  hs’n+1(r/)G(^|r/)dr/ 


1 

I 


r  ext  L 


Ahs’n+1(r ')  -iis’n+1(r') 


^■*4  (0.1.17) 


We  will  separately  analyze  the  integral  term  associated  with  the  exterior  domain, 
Ur,  and  the  integral  term  associated  with  the  impedance  plane,  rext.  The  approaches 
are  similar. 

For  the  first  integral  term,  we  divide  the  exterior  domain  into  a  triangular  mesh 
with  a  sufficient  number  of  triangles  (denoted  M).  Even  though  the  exterior  domain 
will  not  have  to  be  very  large,  because  the  scattered  held  rapidly  diminishes  in  value 
away  from  the  cavity,  it  is  still  important  to  ensure  we  have  a  sufficient  number  of 
triangles  to  ensure  convergence.  Each  triangle,  denoted  Ak,  has  a  centroid,  denoted 
A*,,  and  an  area,  denoted  | | -  We  then  fix  r'  =  Ak  and  approximate  this  integral 
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term  as: 


M 

or  II  us’n+1(r')G(£j\r')dr'  «  a2  Y''  |Afc|  G(^-|Afc)its’n+1(Afe).  (6.1.18) 

/w«  fc=i 


We  note  that  in  approximating  ■us,n+1(A/c)  with  the  Newmark  scheme,  we  have 


us’n+1( Afc)  =  us,n(Ak)  +  5taa-"(Afc)  +  ®  (1  -  2f3)us,n(Ak), 


which  requires  precomputation  of  the  values  of  the  scattered  held  us,n  at  Ak.  Thus, 
using  the  integral  representation  in  (3.2.10),  we  have: 


«s’n(Afc)  = 


=  a 


G{Ak\r')us'n(r')dr' 


UR 


1 

A 

L 


Text  L 


Aua’n+1(r')  -  iis’n+1(r') 


onr, 


(CAMrG4P~^rG-^W. 


dnr / 


dnr> 


Using  the  previous  techniques,  we  can  estimate  us,n(Ak)  as  follows: 


us’n(Ak)^a2  \^n\G(Ak\An)us’n(An)  (6.1.19) 

ahf 

+  a2-us,rt(An)|-ABasex - /  A0 (r)dr  +  |  Afc  |  Grem(Afc  | Afc)  |  (6.1.20) 

v  2  Ztt  ex  J q  J 

+  E 

l=i 

- 1 E  ir“<i  {^«*'”+1(&)  -  «*'n+1te)}  AAM.  (6.1.22) 

The  hrst  term  in  the  estimation  above,  (6.1.19),  is  similar  to  the  previous  exterior 
discretization,  but  it  does  not  include  the  term  when  Ak  =  An.  The  second  term 
(6.1.20)  addresses  this  singularity  by  evaluating  the  hrst  term  of  the  Green’s  function 


k\£j 


dnr, 


,gg( Mi 

dnrt 


(6.1.21) 
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as  previously  described,  and  then  evaluating  the  remaining  terms  of  the  Green’s 
function,  taking  into  account  the  singular  behavior.  In  this  term,  Aease  denotes  the 
length  of  the  base  of  the  triangle,  and  h!  denotes  the  height  of  the  triangle.  The  third 
term,  (6.1.21),  is  the  discretization  along  the  artificial  boundary,  where  un  denotes 
the  value  of  the  scattered  held  as  computed  in  the  interior.  The  fourth  term  (6.1.22) 
is  the  discretization  of  the  impedance  plane. 

Once  us,n( Afc)  is  computed,  we  then  compute: 

h^(Afc)=a2(u^(Afc)-u^(Afc)), 
ua’n(Ak)=tf’n(Ak)  +  6t'yua'n(Ak), 

so  that  -us,n+1(Afc)  and  ■us,n+1(Afc)  can  be  computed  as  desired. 

Similarly,  for  the  second  integral  term  in  (6.1.17),  we  can  divide  Text  into  sufficient 
segments  to  obtain  the  approximation: 

-j[  Aus,n+1(r')  —  iis,n+1(r')  dG^jr  -dx'  « 

y  r ext  l  j  y 

-J  i>2  lrextj  Uus’"+1(&);  -  *s’n+1(6)]  dG&}&.  (6.1.23) 

1=1  l  \  y 

Using  the  same  approach  as  with  the  first  integral  expression,  in  order  to  get 
,~.s;n+i(^)  an(]  ^s,n+ 1(^);  we  again  utilize  (3.2.10)  to  get: 

us’n(&)  «  a2  J]  | An|  G(6|An)us’"(An)  (6.1.24) 

71=  1 

+  E  ISfl, 

3= 1 

-  j  E  lUxU  {***'n+1(a  -  *■-»««„)}  AAM.  (6.1.26) 

m= i  y 

The  hrst  term  in  the  estimation  above,  (6.1.24),  is  similar  to  the  previous  exterior  dis- 


{G(&fe 


dun(^) 

dn„r 


dG{£e\€j 


(6.1.25) 


94 


cretization,  and  accounts  for  all  terms  (where  no  singularities  are  present).  The  second 
term,  (6.1.25),  is  the  discretization  along  the  artificial  boundary,  and  no  singulari¬ 
ties  should  exist  here.  The  third  term  (6.1.26)  is  the  discretization  of  the  impedance 
plane,  and  a  singularity  will  exist  when 

Again,  once  us,n(&)  is  computed,  we  then  compute: 

*',B(e<)  =  a2(«s,B(6)-fi-'B(6)), 

so  that  ■us,n+1(Afc)  and  ■iis,n+1(Afc)  can  be  computed  from  the  prediction  formulas  as 
desired. 

Thus,  the  final  discretized  version  of  the  expression  in  (6.1.17)  is 
* Rus’n+ Htj)  =  S-^a2  |  Afc|  G(£j |  Ak)us'n+1  (Ak) 

k=  1 

-Air-I  U«,’’,+1fe)  -  ftw,+1te)]  aG)fJ6)  }■  (6.1.27) 

71  e=i  L  J  y 

We  note  here  that  S~x  exists  as  a  square,  dense  matrix,  as  previously  described. 
The  two  discretized  integral  expressions  in  the  braces  will  yield  a  Nseg  x  1  matrix, 
where  Nseg  refers  to  the  number  of  segments  used  on  the  artificial  boundary.  When 
multiplied  by  A-1  ,  an  Nseg  x  Nseg  matrix,  this  will  result  in  a  Nseg  x  1  matrix. 
Therefore,  we  multiply  the  result  by  the  “connectivity”  matrix,  [Ih]T ,  (see  (6.1.11)), 
resulting  in  the  desired  iV  x  1  matrix,  where  N  represents  the  number  of  nodes  in 
the  interior  mesh.  This  is  similar  to  Copeland,  Langer,  and  Pusch’s  discretization 
technique  used  in  [20],  and  also  was  referred  to  earlier  in  Jin  in  [55]. 

Hence,  combining  the  discretized  approximations  in  (6.1.12),  (6.1.13),  (6.1.14), 
(6.1.15),  (6.1.16),  and  (6.1.27)  generates  an  approximation  for  the  vector  A”+1  in 
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(6.1.7). 

We  can  concisely  summarize  the  time-stepping  scheme  as  follows  (this  template 
is  adapted  from  Van  and  Wood  in  [104])  : 

1.  Form  the  matrices  K,M,P,  and  Q  as  defined  in  (6.1.3)  ,  (6.1.4),  (6.1.5),  and 

(6.1.6) .  This  constitutes  the  left-hand  side  of  the  matrix  equation  (6.1.2). 

Begin  the  time  loop,  where  n  —  0, 1,  2, . . . ,  N. 

2.  Compute  the  predicted  values  us,n+i  (see  (3.1.3))  and  us’n+1  (see  (3.1.4))  in  the 
interior  Ur.  These  predicted  values  are  based  on  the  finite  element  solution 
computed  in  the  previous  time  step  (see  (6.1.2)). 

3.  Compute  the  predicted  values  us,n+1  (see  (3.1.3))  and  fC,n+1  (see  (3.1.4))  in 
the  exterior  Ur.  These  predicted  values  are  based  on  the  exterior  solution 
computed  from  the  integral  representation  (the  discretized  version  of  (3.2.10)) 
in  the  previous  time  step. 

4.  Construct  the  right-hand  side  vector  F”+1  as  defined  in  (6.1.7). 

5.  Compute  the  solution  of  the  matrix  equation  in  (6.1.2).  This  will  be  the  solution 
in  the  interior  of  the  cavity  Ur. 

6.  Using  the  interior  solution  computed  in  the  previous  step  (un+1),  and  noting 
that  us,n+1  =  un+1  —  umc’n+1  along  Br,  compute  the  solution  in  the  exterior 
domain  from  the  integral  representation  (3.2.10). 

7.  Correct  the  solution  computed  in  the  interior  of  the  cavity  Ur  by  computing 
un+l  (see  (3.1.6))  and  iin+1  (see  (3.1.7)). 

8.  Correct  the  solution  computed  in  the  exterior  Ur  by  computing  us,n+l  (see 

(3.1.6) )  and  us'n+1  (see  (3.1.7)). 
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6.2  Numerical  Results 


We  will  analyze  numerically  two  separate  cavity  geometries:  a  shallow  planar 
cavity,  with  dimensions  1  meter  wide  by  0.25  meters  deep;  and  an  overfilled  cavity, 
with  dimensions  1  meter  wide  by  1  meter  deep,  with  a  semicircle  protrusion  of  0.5 
meter  radius.  Even  though  the  theory  presented  assumes  an  overfilled  cavity,  the 
planar  cavity  results  could  provide  a  good  approximation,  and  more  importantly  will 
provide  a  context  for  the  theory  presented.  It  will  also  be  useful  to  compare  the  electric 
held  depictions  as  well  as  the  RCS  values  for  both  types  of  cavities.  We  previously 
mentioned  that  Wang  provided  numerical  results  for  absorbing  boundary  conditions 
for  curved  boundaries  [108] .  Because  our  assumptions  essentially  approximate  a  first 
order  absorbing  boundary  condition,  this  gives  us  a  baseline  to  help  confirm  the 
accuracy  of  the  depicted  fields. 

We  started  with  the  MATLAB  code  from  Van  and  Wood  in  [104]  as  a  foundation 
for  our  results,  and  adapted  it  to  our  boundary  conditions  and  geometries.  The 
idea  is  to  depict  the  progression  from  the  PEC  case  to  the  IBC  case,  in  which  we 
would  expect  to  observe  the  attenuation  of  the  electric  held  with  a  more  absorbing 
boundary  condition.  For  the  planar  cavity,  we  considered  two  cases:  a  PEC  plane 
(rext)  with  PEC  cavity  walls  (A),  and  a  PEC  plane  with  IBC  cavity  walls.  For  the 
overfilled  cavity,  we  considered  these  two  cases  and  added  a  third  case:  an  IBC  plane 
with  IBC  cavity  walls.  We  also  wanted  to  observe  the  long-term  stability  in  these 
held  depictions,  which  will  validate  the  parameters  chosen  for  the  stability  of  the 
Newmark  scheme. 

We  will  use  the  following  parameters  for  the  Newmark  scheme  for  the  numerical 
simulations:  7  =  0.95,  [3  =  0.5256,  and  St  =  0.0625.  This  ensures  that  Theorem  5.3.1 
is  satisfied  and  the  scheme  is  unconditionally  stable.  We  assume  a  relative  permittivity 
of  er  =  2  for  the  filled  media  in  the  cavity,  which  satisfies  the  assumptions  used  in  the 
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proof  of  Theorem  4.2.3.  We  assume  all  permeability  values  of  y  to  be  1.  We  assume 
the  surface  impedance  values  for  the  IBCs  to  be  rj  =  0.8  if  not  explicitly  mentioned, 
and  we  will  also  use  rj  =  0.2  as  a  comparison  for  a  boundary  condition  approaching  a 
PEC  surface.  Note  that  rj  =  0  represents  a  strict  PEC  surface. 

We  assume  an  incident  Gaussian  pulse  defined  as  follows: 

umc(x,y,t )  =  A—^=  e~r2 , 

1y7[y 

Amplitude 

for 

(t  -  to +  x  cos  6>inc  +  y  sin  0inc) 

r  =  4 - — - . 

T 

In  the  numerical  simulations  for  the  pulse,  we  set  A  =  1,  dinc  =  90°  (=  7t/2),  T  =  2, 
and  to  =  3.  This  time  delay  of  to  =  3  means  at  the  point  (To,  Vo)  =  (0,  0),  the  Gaussian 
pulse  will  reach  its  maximum.  There  are  several  advantages  in  using  the  Gaussian 
pulse,  one  of  which  is  that  it  is  a  good  approximation  of  the  pulse  shape  coming  from 
certain  types  of  lasers  and  other  manmade  systems  [76].  The  depicted  electric  fields 
will  be  the  real  portions  plotted  against  time  as  measured  in  light  meters,  which  is 
the  amount  of  time  for  light  to  travel  one  meter  in  free-space. 

The  Erst  set  of  data  is  for  the  shallow  planar  cavity,  depicted  along  with  the 
incident  wave  and  observation  point  at  (0,0),  in  Figure  8. 

The  first  run  is  depicted  in  Figure  9.  With  the  PEC  plane  and  PEC  cavity  walls  as 
a  benchmark,  we  observe  that  the  case  with  IBC  enforced  at  the  cavity  walls  exhibits 
more  attenuated  characteristics.  We  also  observe  the  stability  of  the  Newmark  scheme 
over  time. 

We  also  want  to  observe  the  effects  as  rj  — *  0,  so  the  subsequent  run  adjusts 
y  =  0.2.  We  would  expect  the  held  to  begin  to  exhibit  the  characteristics  of  a  strict 
PEC  on  the  plane  and  cavity  walls.  This  is  clearly  evident  in  Figure  10  with  the 
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Figure  8.  Shallow  Planar  Cavity  (1  meter  by  0.25  meters) 

second  case  (PEC  plane  and  IBC  cavity  walls)  showing  more  oscillatory  behavior  as 
in  the  strict  PEC  case.  Again,  we  observe  the  stability  of  the  Newmark  scheme  over 
time. 

It  is  essential  to  compute  the  RCS  for  any  scattering  problem,  so  we  generated 
plots  for  two  selected  frequencies  of  289.5  MHz  and  480.45  MHz.  Not  only  are  we 
interested  in  observing  the  effects  of  the  changing  boundary  conditions  on  the  cavity 
models,  but  also  are  interested  in  differences  between  the  planar  and  overfilled  cavities 
as  well.  These  plots  are  monostatic  in  that  the  transmitter  and  receiver  are  collocated 
[8].  They  depict  the  RCS  values  from  a  normal  incidence  angle  of  90  degrees  to  170 
degrees  as  measure  from  the  ground  plane.  We  again  use  Van  and  Wood’s  code  in 
[104]  to  generate  our  results,  in  which  the  time  domain  scattered  held  is  Fast  Fourier 
Transformed  to  obtain  the  frequency  result. 

For  the  Erst  RCS  plot  at  289.5  MHz  in  Figure  11,  we  observe  the  expected  lobing 
pattern  in  both  cases.  We  also  observe  generally  lower  RCS  values  once  the  IBC  is 
introduced,  particularly  for  values  with  an  obervation  angle  less  than  140  degrees. 
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EJV/m) 


Electric  field  at  (0,0)  with  6.  =  90 

'  '  me 


Figure  9.  Shallow  Cavity  -  TM  Solution  at  (0,0)  -  Eta  Point  Eight 
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Electric  field  at  (0,0)  with  6.  =  90 
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Figure  10.  Shallow  Cavity  -  TM  Solution  at  (0,0)  -  Eta  Point  Two 
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TM-RCS  at  289.5  MHz 


Figure  11.  Radar  Cross  Section  at  289.5  MHz 
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In  the  second  RCS  plot  at  480.45  MHz  in  Figure  12,  we  observe  much  more 
apparent  attenuation  of  the  RCS  values  as  the  boundary  condition  is  changed  from 
a  PEC  to  an  IBC,  as  a  significant  gap  exists  between  the  two  cases.  This  higher 
frequency  also  exhibits  the  expected  lobing  patterns,  and  the  PEC  plane  with  IBC 
cavity  walls  exhibits  the  lowest  RCS  values. 

The  second  set  of  data  is  for  the  overfilled  cavity  model  depicted  along  with 
the  incident  wave  and  observation  point  of  (0,  0)  in  Figure  13.  We  point  out  that  the 
selected  observation  point  is  close  to  the  origin  in  order  to  standardize  it  to  the  planar 
case  and  see  if  any  differences  are  observed.  For  these  results,  we  also  include  the  third 
case  of  the  IBC  plane  and  IBC  cavity  walls.  We  also  include  the  mesh  discretization 
used  for  the  interior  domain  (cavity)  (Figure  14)  and  the  exterior  domain  (Figure  15). 

For  the  first  result  in  Figure  16,  we  immediately  notice  the  more  oscillatory  nature 
of  the  electric  field.  This  could  be  clue  to  a  couple  of  reasons:  the  overfilled  cavity 
model  is  physically  deeper  and  has  more  material  surrounding  the  observation  point. 
However,  we  still  observe  the  most  important  result;  that  is,  as  the  boundary  condition 
changes  on  both  the  plane  and  cavity  walls  from  a  PEC  to  an  IBC  we  clearly  see  a 
progressive  attenuation  of  the  depicted  fields.  Also,  this  result  is  truncated  at  50  LM 
for  scaling  and  ease  of  observation.  The  long-term  stability  is  evident  for  the  PEC 
plane  with  IBC  cavity  walls  and  the  IBC  plane  with  IBC  cavity  walls.  For  the  PEC 
plane  with  PEC  cavity  walls,  the  long-term  stability  is  evident  beyond  50  LM  but  is 
not  depicted  here  as  it  is  easier  to  compare  the  three  cases  in  this  view.  The  PEC 
case  is  used  as  a  benchmark  for  comparison,  and  its  stability  was  studied  and  proven 
by  Van  and  Wood  [104], 

As  with  the  planar  case,  we  want  to  observe  the  effects  as  rj  — >  0,  and  this  is 
shown  in  Figure  17.  Even  though  the  differences  between  Figure  16  and  Figure  17 
are  more  subtle,  it  still  shows  that  as  r]  — >  0,  the  fields  approach  the  strict  PEC 
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RCS  (dBm) 


TM-RCS  at  480.45  MHz 


Figure  12.  Radar  Cross  Section  at  480.45  MHz 


104 
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Figure  13.  Overfilled  Cavity  (1  meter  deep;  0.5  meter  radius  of  protrusion) 
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Figure  14.  Interior  Domain  Mesh  -  Overfilled  Cavity 
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Figure  15.  Exterior  Domain  Mesh  -  Overfilled  Cavity 
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EJV/m) 


Electric  field  at  (0,0)  with  Gjnc  =  90 


Figure  16.  Overfilled  Cavity  -  TM  Solution  -  Eta  Point  Eight 


108 


behavior.  This  is  evident  in  that  the  depiction  of  the  held  in  the  PEC  plane  with  IBC 
cavity  walls  increases  in  magnitude  to  follow  the  strict  PEC  case.  It  is  also  evident 
in  the  IBC  plane  with  IBC  cavity  walls  case  in  that  at  approximately  10  LMs,  the 
oscillatory  nature  as  seen  in  the  strict  PEC  case  is  starting  to  form.  The  stability  is 
also  observed  up  to  50  LM,  and  can  be  observed  beyond  50  LM  for  the  PEC  plane 
with  PEC  cavity  walls  case. 

We  also  observe  the  RCS  plots  for  the  overfilled  cavity  model,  keeping  in  mind 
that  we  would  expect  RCS  to  be  affected  due  to  differences  in  shape  and  material. 
For  the  Erst  RCS  plot  at  289.5  MHz  in  Figure  18,  we  note  that  there  is  little  difference 
in  RCS  values  for  all  three  cases  as  they  follow  the  same  pattern.  We  do,  however, 
start  to  observe  the  expected  attenuation  at  observation  angles  beyond  130  degrees. 
We  also  note  the  pattern  is  slightly  different  than  the  planar  case. 

On  the  other  hand,  for  the  RCS  plot  at  480.45  MHz  depicted  in  Figure  19,  we 
see  a  much  clearer  separation  of  the  three  cases.  We  note  the  lobing  and  the  clear 
attenuation  as  the  IBC  is  introduced  on  both  the  plane  and  cavity  walls.  At  this 
higher  frequency,  it  is  also  interesting  to  compare  the  observations  to  the  planar 
cavity  model.  We  note  the  PEC  plane  and  PEC  cavity  wall  RCS  values  are  generally 
lower  for  the  overfilled  cavity  model  than  for  the  planar  cavity  model. 

In  summary,  the  results  of  the  numerical  simulation  clearly  agrees  with  the  the¬ 
ory  presented.  We  observed  the  expected  attenuation  of  the  held  as  the  IBC  was 
introduced  on  the  plane  and  cavity  walls.  We  also  observed  the  changes  in  the  RCS 
values  as  the  IBC  was  introduced.  The  planar  cavity  model  provided  a  good  context 
for  our  theory,  but  ultimately  the  overfilled  model  and  its  corresponding  numerical 
results  indicate  that  our  mathematical  model  can  be  numerically  implemented.  The 
long-term  stability  of  the  models  was  also  verified. 
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Electric  field  at  (0,0)  with  0.  =  90 
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Figure  17.  Overfilled  Cavity  -  TM  Solution  -  Eta  Point  Two 
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TM-RCS  at  289.5  MHz 


Figure  18.  Overfilled  Radar  Cross  Section  at  289.5  MHz 
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TM-RCS  at  480.45  MHz 


Figure  19.  Overfilled  Radar  Cross  Section  at  480.45  MHz 
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VII.  Conclusions  and  Future  Work 


7.1  Conclusions 

In  conclusion,  given  an  incident  electromagnetic  wave  impinging  on  an  overfilled 
cavity  in  an  impedance  ground  plane,  we  have  determined  the  resulting  scattered 
fields.  Moreover,  we  have  proven  that  this  problem  is  well-posed  and  can  be  im¬ 
plemented  through  the  use  of  the  FEM.  We  have  established  a  valid  mathematical 
model  for  detecting  scattered  waves  from  cavities.  In  particular,  we  derived  an  inte¬ 
gral  representation  and  Green’s  function  for  the  solution,  and  used  integral  operator 
theory  in  establishing  well-posedness.  A  key  step  in  the  uniqueness  and  existence 
proofs  was  analyzing  the  properties  of  the  Steklov- Poincare  operator.  This  was  im¬ 
portant  because  an  exact  series  solution  does  not  exist  for  IBC  conditions,  and  using 
already  established  properties  of  integral  operators  helped  facilitate  this  proof.  One 
advantage  of  this  approach  was  that  we  avoided  the  use  of  the  hypersingular  opera¬ 
tor,  which  can  be  difficult  to  numerically  evaluate,  particularly  at  the  singularities. 
Nevertheless,  we  noticed  an  elegant  connection:  the  mathematical  formulation  sets 
the  foundation  for  the  variational  formulation.  This,  in  turn,  is  the  basis  for  the  finite 
element  analysis,  which  is  used  for  the  numerical  simulation.  We  observed  that  not 
only  changing  boundary  conditions  but  changing  cavity  geometries  had  an  effect  on 
the  scattered  field  and  the  RCS  values. 

We  also  have  two  additional  contributions.  In  constructing  the  finite  element 
approximation,  we  provided  a  method  for  numerical  implementation  that  avoids  use 
of  the  hypersingular  operator.  This  approach  in  the  literature  seems  limited  and 
almost  nonexistent.  Though  the  numerical  simulation  involved  a  direct  analysis  and 
evaluation  of  the  equations  of  the  integral  representation,  this  method  may  be  more 
efficient  when  fully  implemented.  Also,  we  were  able  to  apply  and  refine  previous 
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research  for  parameter  requirements  to  ensure  unconditional  stability  for  arbitrary 
time  steps  for  our  application  of  the  Newmark  scheme. 

One  of  the  most  important  observations  involved  the  discretization  of  the  Newton 
potential  terms,  one  of  which  was  the  infinite  space,  Ur,  above  the  cavity  opening 
(see  (3.4.8)).  We  previously  mentioned  that  this  term  results  from  the  inhomogeneity 
in  the  PDE,  and  just  requires  proper  discretization,  as  the  finite  element  method  is 
not  used  in  this  portion.  During  the  numerical  implementation,  Ur  was  relatively 
small  (2  meters  wide  by  1  meter  long;  see  Figure  15),  as  the  scattered  held  becomes 
quite  small  and  negligible  outside  of  this  region.  However,  the  real  breakthrough  in 
the  research  came  when  it  was  discovered  that  the  solution  did  not  converge  until  a 
sufficient  number  of  elements  were  generated  for  this  region.  We  believe  this  might  be 
due  to  the  fact  that  we  were  working  with  the  constant  a2,  which  had  an  extremely 
large  value  of  almost  500.  As  a  result,  when  discretizing  this  portion  of  the  Newton 
potential,  even  though  the  finite  element  method  or  adaptive  mesh  refinements  are  not 
used,  this  refinement  may  have  helped  scale  the  problem  and  avoid  solution  “blow  up” 
at  the  appropriate  steps.  This  is  an  important  consideration  in  future  applications 
and  extensions  of  this  work  in  the  time  domain. 

7.2  Future  Work 

It  is  obvious  that  this  research  has  a  lot  of  applications  and  opens  up  many  avenues 
for  future  work.  Considering  our  general  problem  statement,  we  only  investigated  the 
TM  polarization  case,  but  the  TE  polarization  case  can  be  explored  as  well.  Also,  this 
problem  can  be  translated  into  three  dimensions,  which  would  not  only  increase  its 
applicability,  but  also  be  easier  to  work  with  in  terms  of  a  simpler  Green’s  function. 
This  was  a  direct  scattering  problem,  but  the  IED  detection  application  lends  itself 
well  to  an  inverse  scattering  problem.  That  is,  given  a  scattered  wave,  determine  the 
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properties  of  the  scattering  object.  Inverse  scattering  theory  is  a  popular  and  growing 
area  of  research. 

Considering  the  numerical  analysis  of  the  problem,  we  only  generated  the  real 
portion  of  the  scattered  held,  but  future  research  could  also  look  at  the  imaginary 
portions  of  the  held  generated  due  to  the  nature  of  the  Green’s  function.  Also, 
RCS  plots  could  be  analyzed  at  several  other  frequencies.  Several  of  the  parameters 
could  also  be  modified  and  analyzed  for  their  effects  to  include:  different  values 
of  er  (to  include  real  and  imaginary)  to  rehect  different  types  of  cavity  material; 
different  values  of  the  surface  impedance  77;  changing  the  angle  of  incidence  to  see 
its  effects  on  the  scattered  held;  changing  the  observation  points  along  and  within 
the  artihcial  boundary;  different  overfilled  cavity  geometries  (shallow  versus  deep 
cavities);  changing  the  stability  parameters  used  in  the  Newmark  scheme  to  include 
differing  time  steps;  and  considering  the  effects  of  a  Gaussian  pulse  versus  a  continuous 
wave.  Obviously,  there  are  many  areas  here  that  could  be  expanded  and  analyzed. 

We  mentioned  that  the  IBC  assumption  we  developed  approximated  a  hrst  order 
boundary  condition,  but  a  more  robust  assumption  could  be  made  involving  the  use 
of  a  Fourier  transform.  The  Green’s  function  itself  for  an  IBC  can  be  explored  further, 
as  Hein  did  thoroughly  in  his  dissertation  [42],  It  was  clear  from  the  literature  that 
there  is  little  consistency  in  the  Green’s  function  expressions  for  an  IBC.  There  are 
multiple  representations  for  different  terms,  and  each  are  quite  detailed.  The  main 
difficulty  seems  to  be  in  numerical  implementation  and  the  evaluation  of  singularities 
when  using  a  boundary  integral  method.  More  standardization  in  this  area  would  be 
very  beneficial. 

As  mentioned  during  the  finite  element  analysis,  we  only  provided  a  priori  esti¬ 
mates,  but  adaptive  mesh  refinement  could  be  used  to  generate  a  posteriori  estimates, 
or  higher  polynomial  degree  basis  functions  could  even  be  used  to  help  refine  the  so- 
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lution.  For  the  finite  element  approximation  we  describe,  the  idea  was  to  use  it  as 
a  foundation  for  the  numerical  implementation.  However,  we  built  upon  an  existing 
code  which  more  directly  approximated  the  integral  representations.  The  advantage 
was  that  this  did  not  require  inverse  matrix  operations,  but  the  disadvantage  was  that 
the  calculations  involved  hypersingular  expressions  with  multiple  derivatives,  where 
singularities  were  difficult  to  assess.  The  finite  element  approximation  we  describe 
could  be  implemented  as  it  avoids  use  of  hypersingular  expressions,  and  could  possibly 
be  more  stable  and  accurate. 

Finally,  it  is  without  question  that  this  research  has  opened  several  areas  for 
further  exploration  that  could  greatly  benefit  not  only  the  academic  discipline  but 
also  the  military.  Given  the  importance  of  radar  applications  and  the  importance  of 
object  detection  and  analysis  in  today’s  military,  this  research  can  be  applied  in  many 
ways.  Clearly,  then,  this  work  sets  the  foundation  for  extremely  important  topics  to 
today’s  warfighter  on  the  battlespace. 
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Appendix  A.  Derivations 


Because  we  assume  a  constant  surface  impedance  for  our  impedance  boundary 
conditions  in  the  time  domain,  this  assumption  resembles  a  first-order  absorbing 
boundary  condition  according  to  Jin  (see  [55])  : 
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14.  ABSTRACT 

We  consider  the  transient,  or  time-domain,  scattering  problem  of  a  two-dimensional  overfilled  cavity  embedded  in  an 
impedance  ground  plane.  This  problem  is  a  significant  advancement  from  previous  work  where  more  simplified  boundary 
conditions  were  used,  which  can  limit  the  number  of  applications.  This  research  supports  a  wide  range  of  military 
applications  such  as  the  study  of  cavity-like  structures  on  aircraft  and  vehicles.  More  importantly,  this  research  helps 
detect  the  biggest  threat  on  today’s  battlefield:  improvised  explosive  devices.  An  important  step  in  solving  the  problem 
is  introducing  an  artificial  boundary  condition  on  a  semicircle  enclosing  the  cavity;  this  couples  the  fields  from  the  infinite 
exterior  domain  to  those  Helds  inside.  The  problem  is  first  discretized  in  time  using  the  Newnrark  scheme,  and  at  each 
time  step,  we  derive  the  variational  formulation  and  establish  well-posedness  of  the  problem.  This  sets  the  foundation  for 
the  finite  element  method  used  in  the  numerical  analysis.  Using  both  a  planar  and  overfilled  cavity  model,  we  provide 
numerical  results  through  the  depictions  of  the  scattered  electric  field  and  radar  cross  section  of  the  cavities. 
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